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APPLICATION OF SYSTEM IDENTIFICATION TO ANALYTIC ROTOR MODELING 
FROM SIMULATED AND WIND TUNNEL DYNAMIC TEST DATA 

Part II of Final Report under Contract NAS2-7613 


Abstract 

This report begins with an introduction to aircraft state and parameter 
identification methods. A simplified form of the Maximum Likelihood method 

• i 

is selected to extract analytical aeroelastic rotor models from simulated 
and dynamic wind t unn el test results for accelerated cyclic pitch stirring 
excitation. The goal is to determine the dynamic inflow characteristics 
for forward flight conditions from the blade flapping responses without 
direct inflow measurements. The rotor blades are essentially rigid for 
inplane bending and for torsion within the frequency range of study, but 
flexible in out-of-plane bending. Reverse flow effects are considered for 
high rotor advance ratios. 

Two inflow models are studied; the first is based on an equivalent 
blade Lock number, the second is based on a time delayed momentum inflow. 

In addition to the inflow parameters, basic rotor parameters like the blade 
natural frequency and the actual blade Lock number are identified together 
with measurement bias values. The effect of the theoretical dynamic inflow 
on the rotor eigenvalues is studied. A relation between the accuracy of 
the identified parameters and the length of the input data is established in 
simulation studies. 

It is found that the first inflow model using an optimized equivalent 
blade Lock number is very accurate for rotor advance ratios of .4 and above, 
while for lower advance ratios, the second inflow model using a time 
delayed momentum inflow provides better accuracy. For the first inflow 
model the identified equivalent Lock number deviates systematically from 
theoretical values established in the literature. The identified analytical 
models are verified by predicting the test results not used in the 
identification process. 



Preface to Pinal Report under Contract NAS 2- 76 13 


Work under Contract NAS2-7613 started on July 1, 1973. The contract 
was originally awarded for a 3 year period. 

Due to the slower than anticipated progress of the experimental work, 
not all research goals had been achieved by 30 June 1976. Since less than 
the anticipated cost for personnel and equipment had been spent, the 
research contract was extended by a year without increase in funding. 

The research goals as stated in the contract were: 

(a) Assess analytically the effects of fuselage motions on stability 
and random response. The problem is to develop an adequate but not 
overly complex flight dynamics analytical model and to study the 
effects of structural and electronic feedback, particularly for 
hingeless rotors. 

(b) Study by computer and hardware experiments the feasibility of ade- 
quate perturbation models from non-linear tri'i conditions. The 
problem is to extract an adequate linear perturbation model for the 
purpose of stability and random motion studies. The extraction is 
tc be performed on the basis of transient responses obtained either 
by computed time histories or by model tests. 

(c) Extend the experimental methods to assess rotor wake-blade 
interactions by using a 4-bladed rotor model with the capability 
of progressing and regressing blsde pitch excitation (cyclic pitch 
stirring), by using a 4-bladed rotor model with hub tilt stirring, 
and by testing rotor models' in sinusoidal up or side flow. 
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Including the final report, 10 reports under Contract NAS2-7613 
have been submitted. They are listed as P. 1 to P. 10 at the end of 
the Preface. P. 1 and P. 10 pertain to research goal (a). P. 2, P. 4, 

P. 6, P. 7, P. 8, P. 9, pertain to research goal (b). P. 3 and P. 5 
pertain to research goal (c). The latter is not as yet complete since 
neither hub tilt stirring nor testing is sinusoidal up or side flow 
has been performed. While P. 10 describes only work done during FT 1977, 

P. 8 and P. 9 combine both FT 1977 work results and summaries of earlier 
results, so that the three parts of the Final Report can be read without 
recourse to the earlier reports. P. 8 includes much new material not 
available when the preceding Tearly Report P. 7 was written. The 
experimental data of P. 9 have all been obtained in FT 77. 

So far 3 publications came out of the research under Contract NAS2-7613. 
They are listed as P. 11, P. 12, P. 13. 

Lxst of Reports and Papers 
under Contract NAS2-7613 

P 1. Hohenemser, K. H. and Tin, S. K. , "Methods Studies Toward Simplified 
Rotor-Body Dynamics", Part I of First Tearly Report under Contract 
MAS2-7613, June 1974. 

P 2. Hohenemser, K. H. and Tin, S. K. , "Computer Experiments in Pre- 
paration of System Identification from Transient Rotor Model 
Tests", Part II of First Tearly Report under Contract NAS2-7613, 

June 1974. 

P 3. Hohenemser, K. H. and Crews, S. T. , "Experiments with a Four-Bladed 
Cyclic Pitch Stirring Model Rotor", Part III of First Tearly Report 
under Contract NAS2-7613. 

P 4. Hohenemser, K. H. , Banerjee, D. and Tin, S. K., "Methods Studies 
on System Identification from Transient Rotor Tests", Part I of 
Second Tearly Report under Contract NAS2-7613, June 1975. 
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P 5. Hohenemser, K. R. and Crews, S. T. , "Additional Experiments with a 
Four-Bladed Cyclic Pitch Stirring Model Rotor". Part II of 
Second Yearly Report under Contract NAS2-7613, June 1975. 

P 6. Hohenemser, R. H., Banerjee, D. and Yin, S. K. , "Rotor Dynamic 
State and Parameter Identification from Simulated Forward Flight 
Transients", Part I of Third Yearly Report under Contract 
NAS2-7613, June 1976. 

P 7. Hohenemser, R. H. and Crews, S. T., "Rotor Dynamic State and 

Parameter Identification from Hovering Transience" , Part II of 
Third Yearly Report under Contract NAS2-7613, June 1976. 

P 8. Hohenemser, R. H. and Crews, S. T. , "Unsteady Hovering Rotor 

Wake Parameters Identified from Dynamic Model Tests", Part I of 
Final Report under Contract NAS2-7613, June 1977. 
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Contract HAS2-7613, June 1977. 

P 10. Hohenemser, K. H. and Yin, S. K. , "Finite Element Stability 
Analysis for Coupled Rotor and Support Systems", Part III of 
Final Report under Contract NAS2-7613, June 1977. 

P 11. Hohenemser, K. H. and Yin, S. K. , "On the Use of First Order 

Rotor Dynamics in Multiblade Coordinates", 30th Annual National 
Forum of the American Helicopter Society, May 1974, Preprint 831. 
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APPLICATION OF SYSTEM IDENTIFICATION TO ANALYTIC ROTOR 
MODELING FROM SIMULATED AND WIND TUNNEL DYNAMIC TEST DATA 


1. INTRODUCTION 


System Identification is a method of correlating a mathematical 
model of a system with transient responses obtained either 
experimentally or from the time history of a more complete analytical 
model of the system. It is particularly useful if a linear 
perturbation model of a basically non-linear system is to be identified. 

Methods for state and parameter estimation from transients are 
widely used in aircraft testing (1)*, (2), (3) and (4). The problem 
is to obtain optimum estimates (based on certain performance criteria) 
of initial states and of unknown parameters (derivatives) from noisy 
measurements of some inputs and response variables. In most cases 
of airplane parameter identification a constant coefficient system 
is used as an analytical model. For lifting rotor applications, 
a periodic coefficient system model may be required (5) . 


*The numbers in parantheses in the text indicate references in 
the Bibliography. 
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While the identification of stability and control derivatives for 
fixed wing and rotary wing aircraft from transients is by now a 
well established method, the question arises whether or not similar 
methods can also be used to obtain some insight with respect to 
aeroelastic rotor characteristics. When using linear perturbation 
equations from a non-linear trim condition, frequency response 
testing is (me way to correlate a mathematical model with 
experimental dynamic data. This method has been occasionally used 
for wind tunnel rotor model tests (6) and (7) . A less laborious 
and less time consuming method of rotor dynamic testing is to 
extract analytical perturbation models from transient rotor responses. 
The study described herein is the first attempt of accomplishing 
this objective for an aeroelastic rotor model in forward flight 
conditions. Since new ground is covered by extending aircraft 
identification methods to aeroelastic rotor problems, extended 
simulation studies have to be performed to assure the feasibility 
of the identification process. 

A lifting rotor is both structurally and aerodynamical ly a 
highly complex system that has not as yet been fully explored. 

Each blade has natural modes with predominantly out-of-plane 
(flapping), inplane (lead- lag) and torsional motions that are 
structurally and aerodynamical ly coupled. Furthermore the various 
blades of a rotor are also coupled by hub angular or linear 
motions, by control element motions and by the rotor inflow. 

In the following studies a drastically simplified analytical rotor 
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model will be used, where the blades are essentially rigid in inplane 
and torsional motion and where the rotor hub is also rigidly 
supported. Thus only out-of-plane (flapping) blade motions are 
considered and the only inter-blade coupling is from the rotor 
inflow. While there is considerable literature on blade flap- 
bending. for example (8) , (9) and (10) , and while the steady 
rotor inflow has been frequently studied, for example in (9) , 
the only dynamic inflow theory applicable to forward flight conditions 
is given in (11) . The present study uses for moderate rotor 
advance ratio an analytic inflow model that is an extension of that 
given in (11). It correlates this model with transient wind 
tunnel test data with the help of state and parameter identification. 
The corresponding study for hovering conditions is presented in (12) . 

In addition to the inflow model of (11) a substantially simpler 
inflow model will also be studied, based on the replacement of the 
blade Lock number by an optimized equivalent value. This concept 
was originally suggested for steady rotor conditions in (8) and for 
dynamic rotor conditions in (13). Theoretically the equivalent Lock 
number should be a complex number but will be assumed here as a real 
number, corresponding to a quasi-steady analysis. The results with 
the two selected inflow models will be compared to the rotor 
responses when the dynamic inflow is entirely neglected as is done 
in most current aeroelastic rotor analyses. 
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It is seen in (14) that in hovering conditions and using the 
theoretical inflow model of (11), the damping of the regressing 
rotor flapping mode can be substantially reduced, particularly 
at low collective pitch setting. The effect on rotor eigenvalues, 
of the inflow models studied herein for forward flight conditions, 
will be determined to find out the applicability of the various 
inflow models and the frequency ranges in which they are suitable. 
The sensitivity of the rotor eigenvalues to variations in the 
parameters will also be considered. 
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2. AIRCRAFT STATE AND PARAMETER IDENTIFICATION METHODS 

The review of identification methods to be given in the 
following is by no means complete. Only the most important methods 
are discussed. Only rough outlines are given for the various 
methods. Details of the derivations and of the application 
algorithms are found in the cited literature. 

2.1 ELEMENTS OF SYSTEM IDENTIFICATION FROM TRANSIENTS 

System identification is the process of extracting numerical 
values for system parameters and other subsidiary parameters 
(process and measurement noise covariances, bias, initial states, 
etc.) from the time history of control or other inputs and of the 
resulting system responses. A schematic for the measurements is 
shown in Figure 1. The process of system identification involves 
five steps: 

1. Selection of a suitable input that insures participation of 
all important modes of the system in the transient response. 

2. Selection of sufficiently complete and accurate instrumentation 
to measure the key input and output variables. 

3. Selection of a mathematical model that adequately represents 
the actual system characteristics. 

4. Selection of an efficient criterion function and estimation 
algorithm for the identification of the unknown system 


parameters. 
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EXTERNAL DISTURBANCE 



Figure 1. Schematic of Measurements for System Identification* 
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Figure 2. Illustration of System Identification. 
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5. Validation of the identified mathematical model by comparing 

its results to test results not used for the system identification. 
The concept of system identification is illustrated in Fxgure 2. 
The design input is fed both to the actual system and to its mathe- 
matical model that contains the unknown parameters. The measured 
response, polluted by measurement noise, is compared with the 
computed response from the mathematical model. The difference 
between these two responses, the response error, is used in the 
parameter estimation technique based on the criterion function and 
optimizing technique. The estimation algorithm may also use 
apriori information, e.g., initial statistics of the parameters. 

Here we will be mainly concerned with the fourth of the previously 
listed steps, that is with the various estimation algorithms. 

The mathematical representation of the system will be given 
in the non-linear case by: 

System equation x(t) = f(x,u,t) + r(t)w(t) ( 1 ) 

Initial condition x(t = 0) = x Q 

Measurement Equation y(t) = h(x,u,t) + v(t) , 2 ) 

If the system is linear, equations 1 and 2 reduce to 
x(t) = F(t) x (t) + G(t) u(t) + r(t) w(t) 

y(t) = H(t) x (t) + D(t) uft) + g+ v(t) '^) 
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2.2 CLASSIFICATION OF IDENTIFICATION ALGORITHMS 

The various estimation algorithms can be classified into 
two groups presented in Table 1. The first group listed in Table 1 
above the double line is based on statistical regression and does 
not admit a probabilistic interpretation. The algorithms listed 
in Table 1 below the double line are based on probabilistic 
interpretation. In the equation error estimate no measurement 
noise is modeled; the following 4 methods include both measure- 
ment and system noise, while in the output error estimate no 
system noise is modeled. The various algorithms listed in Table 
1 will be discussed in the following sections. 

2.3 EQUATION ERROR ESTIMATES 

Fquation error methods assume a performance criterion that 
minimizes the square of the equation error (process noise) . 

They are least squares techniques and they require the knowledge 
of all response variables (states) and their derivatives. In 
the so called least squares method the unknown parameters are 
selected such that the integral over the square of the state 
equation error is minimized, see for example (2) . With 
equation 1 we have the error function (the upper integration limit 
T is the time over which the measurements are taken) 


* 

/ 


fx - f(x,u,e,t)] r w[x-f (x,u,e,t)]dt 


j 


(5) 



Table 1. Classification of Estimation Algorithms. 
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V is a positive definite weighting matrix. An appropriate 
choice for W would be Q" * where Q is the covariance of the 
process noise. For the usual digital data processing, the 
variables x, x, u are sampled, and only available at discrete 
time points t^. Mathematically the sampling process can be 
expressed by multiplying the system equation with the delta 
function 6(t - t^). The integral of equation 5 then becomes 
a sum. One can use instead of the delta function also a different 
"method function", for example e" st , that would allow taking the 
Laplace transforms. 

If the system is linear in the unknown parameters 9, the 

system equation can be written in the form 

x(t ) * F(x,u,t>8 ♦ T(t) w(t ) (6) 

and the performance criterion (5) becomes 
T 

J a f tx - F(x,u,t)0] T W[x - F(x,u,t)G]dt (7) 

Since the function inside the integral has continuous deriva- 
tives with respect to 0 we set 

3J/30 = 0 (8) 


thus resulting in the closed form sclution 


T - T 

£ J F T (x,u,t)W FU.u.tJdtJ ' 1 J F T (x,u,t)H x(t)dt 


( 9 ) 



11 


The first factor is the covariance matrix of the estimate. 

If the system is non-linear in the unknown parameters, the solution 
equation 9 can be replaced by an iterative solution where F(x,u,t) 
is substituted by dfCXiU.O^tVdQ and § on the left hand side is 

A 

replaced by . It can be shown (for example (1)) that the 

parameters in the n" 1 row of f(x,u,0,t) are independent of all 
the elements of x(t) except x n (t). This independence is one of the 
drawbacks of the least squares method, in that only one of the 
measured state derivatives is used in determining a given row of 
the f(x,u,9,t3 matrix. If one of the signals has not been measured, 
the least squares method does not provide an estimate of the 
parameters related to that signal. This independence also illustrates 
the fact that the estimate of one row of the f(x,u,Q,t) matrix 
is obtained independent of the other rows, and no "trade-off" 
can be wade between elements in different rows to improve the 
estimate. 

For some applications it is practical to include the state 
vectors in the error minimization. In the modified least squares 
method a combination of the standard least squares with the 
integrated least squares is used. The parameters obtained by this 
method not only trace the derivative of the state but also the 
state itself over the selected time interval. The performance 
criterion now includes in addition to the equation error also the 
integrated equation error: 
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J * J | |x(t ) - F(x, 


t 

u,t)0 ♦ J X 


(r)dT - 


t 

J F(x,u,T)e dt|| w dt (10) 


where W is a positive definite weighting matrix and where 

,T 


I |a| | w a a 1 w a 


(ID 


Minimizing the expression, equation 10 results in the estimate 
T t 


J (F(x,u,t) + J F(x,u 


,t) dt} W{F (x,u,T ) + 


t 

£ F(x,u,x) dt} 


dt 


-1 


( 12 ) 


[/ 


T t t 

(F(x,u,t) + J* F (x,u, r )dT}^W{x ( t ) + J x ( t ) dt)dt 


This method has the same row independence off (x,u,0,t) as the 
standard least squares method. 

Since these methods do not allow for measurement errors, they 
result in biased estimates when this type of error does exist. 

When measurement errors are small, as is increasingly the case 
in modem instrumentation, the equation error method *is preferable 
over other methods because of its simplicity. It is widely used 



13 


also when measurement errors are substantial and then serves as 
start-up technique for the output error and other iterative 
methods . 

In many applications, measurements of some of the responses 
or their derivatives are not available. If the response but not 
the rate of response is measured, it is tempting to differentiate 
the measured response. However, the differentiation of measured 
data introduces additional uncertainty so that this technique is 
usually inaccurate. If e is used as a methods function, 

Laplace transforms can be used. The estimation then reduces to 
an algebraic manipulation of the data that avoids their differen- 
tation. The Laplace transform technique as a substitute of 
differentiating measurement data is discussed in (16) . 

2.4 BAYESIAN AND QUASI -BAYES IAN PARAMETER ESTIMATES 

In the preceding methods we specified a cost criterion J 
that represented the "loss" resulting from an incorrect estimation 
of the unknown system parameters. The parameters were then 
selected in such a way as to minimize the loss. If a priori 
probabilities exist not only for the measurement errors but also 
for the unknown parameter vector 9 then one can define an expected 
loss and select the parameter vector in such a way as to minimize 
this expected loss. Such an estimate is called a Bayesian estimate, 
see for example (17). 

The form of a Bayesian estimate depends on the form of both 
the loss function and of the a priori probability distribution of 
the measurement and the parameter vector. 
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For the particular case of positive semi -definite quadratic loss 
functions, the Bayesian estimate is the mean of 6 conditioned on 
the observations. This is true regardless of the distribution of 
measurement and parameter vector (18) and (19). It has also been 
shown that for the case of unimodal symmetric a posteriori distri- 
bution of the parameters given the observations, the Bayesian 
estimate is the conditional mean for all loss functions which 
are symmetric and convex upwards. For these reasons the Bayesian 
estimate can be defined generally as the conditional mean of the 
parameter distribution. 

In order to compute the conditional mean, it is first necessary 
to determine the conditional probability density for 9. This 
density can be written from Bayes rule as (Z is the set of all 
observations) 

P (9/Z) = p (2/0) p (0)/p(Z) (13) 

The denominator is a normalizing factor determined from 

P (3) = J P(Z/9) p(0) de (14) 

all 0 

The optimal Bayesian estimate is now given by 

e = (i/p(z) ) /e p (z/ 0 ) p(e) de (is) 

all e 

In general, the evaluation of equations 14 and 15 would require the 
solution of the system equations for all possible values of the 
parameter vector 0. 
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This is a large effort, especially if the dimension of 0 is large. 

If p(0/Z) is unimodal and symmetric about its mean value, the 
conditional mean corresponds to the mode. Since p(Z) is merely 
a scale factor the finding of the mode requires neither the 
evaluation of the integral in equation 14 nor that in equation 15. 

The mode 0 of 0 has the property 

max p (8®) - max p(Z/8) p(0) = p(Z/0) p(0) (16) 

0 

Even if the a priori density p(0) is symmetric it does not follow 
that the conditional density p(6/Z) is also symmetric since in 
general the observations depend non- linearly on the parameters. 
Estimation according to equation 16 is, therefore, called 
"quasi-Bayesian" estimation. Another designation used for example 
in (3) is maximum a posteriori probability (MAP) parameter estimate. 
Since the logarithm is a monotonic function of its argument, we can 
replace equation 16 by maximizing the expression 

J = log P(Z/0) + log P(0) (17) 

If no a priori information about the parameters 0 is available, 
that is, if the a priori density is uniform, p(9) = constant, the 
quasi-Bayesian estimate reduces to the "Maximum Likelihood" estimate 
which involves finding the maximum of p(Z/9) . 

2.5 ESTIMATES ASSUMING GAUSSIAN DISTRIBUTIONS 

The evaluation of equation 17 becomes particularly convenient 
if we assume Gaussian densities for the parameters, for the observations 
and for the system states. In linear systems and linear measurement 
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equations (equations 3 and 4) one needs only to assume that the 
system noise w(t) and the measurement noise v(t) is Gaussian. 

It then follows that states x(t) and observations y(t) are also 
Gaussian. For non-linear systems with Gaussian noise, p(Z/0) 
tends to a Gaussian density as the sampling rate is increased (see 
for example (1), p.29). The assumption of Gaussian densities for 
all variables is, therefore, a reasonable one. Since 0 is a m x 1 
vector we now have the a priori density 

f(0) * | p 0 1 ” 1/2 (2wr m/2 exp - (i/2)(e-e) T p" 1 (e- 0 ) (is) 

Except for a constant additive term, log p(0) is now given by 

logp(e) * (-1/2) (0-0 ) T P" 1 (0-0) (19) 

In order to obtain an expression for log p(Z/9) in equation 17, 
we assume that Z consists cf N consecutive observations y ( 1) .. y(N) . 

Z * Y n = (y(l), . . . y(N)> (20) 

With successive application of Bayes rule we obtain 

P(Y N /0) = P(y(l), . . y(N)/0) = p (yOO/Y^,©) p(Y, J _ 1 /0) 

N 

= it p(y(j)/v. .,0) (21) 

j=i 

Taking the logarithm we have 

N 

log p (Y M /0 ) = E log p(y( j )/Y . , 0) (22) 

i=i 3 

( y(j)/Yj j,9) is the observation estimate at time j given all 
preceding observations and given the parameters. We denote the 
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observations by y(j) and its expected value and covariance 

A 

respectively by y(j/j-l) and B(j/j-l). We further denote the 
"innovation" by 

y(j ) - y( j/j-1) * v(j) (23) 

Since y(j) is a r x 1 observation vector, its Gaussian 
density is 

p(y( j )) = |B(j/3-l)|- 1/2 (2w)“ r/2 exp|-(l/2)v T (j) 

(24) 

B' 1 (j/j-l)v(j)| 

Taking the logarithm of equation 24, summing according to 
equation 22, inserting in equation 17 and inverting the sign we have 
now to minimize the expression (see also equation 19) 

M T 1 

E (v (j)B A (j/j-l)v(j ) + log |B( j/j-1) | } 

1=1 

+ o-e) T p" 1 (e-¥) (25) 

If no a priori information is available before taking observations, 
the last term in the expression 25 is constant and we then have 
the criterion for the Maximum Likelihood estimation. Bayesian or 
quasi-Bayesian estimation is rarely used since a priori densities 
for the parameters are in most applications not available. 

2.6 MAXIMUM LIKELIHOOD PARAMETER ESTIMATION 

According to expression 25, Maximum Likelihood estimation is 
equivalent to minimizing the so-called likelihood function 

N T , 

= Z [v (j)B A (j/i-l)v(j) + log | B( j / j-1 ) I ] 

1=1 


J(0) 


(26) 
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In the presence of system noise the mir.imicaticn of the 
expression 26 is very difficult. When going from time j-1 to 
time j one first has to solve the prediction equations for the 
estimate of the state and for its covariance. Assuming the linear 
system equation 3 with zero mean Gaussian system noise w(t) the 
prediction is given by 

£(1/1-1) = F *0/1-1) + 1 u(t) , (1-1) i t < j (27) 
p'1/j-i) = r P(j/i-i) ♦ p(1/i-df t + r q r T (28) 

where Q is the system noise covariance and P the state covariance. 
These equations use the estimated state and its covariance at time 
j-1: x(j-l/j-l) and P(j-l/j-l), to predict the state and its 

covariance at time j: x(j/j-l) and P(j/j-l). This is the pre- 

diction before we know the result of the observations at time j. After 
the observations y(j) have been made the optimum estimate is given 
by the Kalman filter equations for the state and for its covariance: 

*(1/D = *<1/j-l) * K j) (yd) - H *0/1-1)) (29) 

P(1/i) = (I - K ( 1 ) H) P(1/j-l) (30) 

with the filter eain 


K(1 ) - P ( i / 1 - 1 ) H r (H P(j/j-l) H T + R)” 1 


(31) 
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The covariance of the observations B(j/j-l) that occurs in the 
cost function 26 is given in terms of the state covariance before 
observations by 

B(j/j-l) * H P(j/1-1) H T + R (32) 

Thus the terms in the expression 26 that is to be minimized require, 
the solution of the prediction equations 27 and 28 for each time 
interval and of the up-date equations 29 and 30 at each sampling 
time together with the solution of the measurement equation 4. 

(1) gives an algorithm for the solution of the problem. However, 
due to its complexity this algorithm has not as yet been applied 
to a practical problem of aircraft parameter estimation, see for 
example (20) . 

The problem of minimizing the expression 26 is greatly simpli- 
fied if the observation covariance B(j/j-l) can be assumed 
constant. This is for example true for zero system noise, when 
according to equation 32, B(j/j-l) = R. The problem then reduces 
to minimizing the cost function 

H T -1 

J(e) » l v (j) R v(j) (33) 

j=l 

where v(j) is given by the innovation term 23. Since equation 33 
represents (according to equations 2 and 4) the sum of the measure- 
ment error squares, the estimation with equation 33 is also called 
output error method of estimation. There are several algorithms 
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available to perform the optimization of J(0) from equation 33. 

The most widely used is the modified Newton-Raphson or quasilinear- 
ization method. It has the advantage that the sensitivity or information 
matrix is obtained as a byproduct. The inverted information matrix 
gives the Cramer- Rao lower bound for the parameter covariance. 

This lower bound is found in many applications to be a more meaningful 
measure of the accuracy of the parameter estimate than the parameter 
covariance obtained from the equation error method (first factor in 
equation 9) . 

A block diagram of the complete and the simplified Maximum 
Likelihood identification procedure is shown in Figure 3. The iteration 
loop is indicated by double lines. Neglecting the three signals shown 
by dashed lines, the Kalman filter reverts to the deterministic 
solution of the system equations. 

2.7 SOME PROPERTIES OF THE MAXIMUM LIKELIHOOD ESTIMATE 

The Maximum Likelihood estimation technique has several 
theoretically justifiable properties which makes it the best accepted 
estimation technique to date. Some of the proven properties of the 
Maximum Likelihood method are: 

1. The Maximum Likelihood estimate is consistent, i.e., the parameter 
estimates converge (in probability) to its true values as the number 
of observations N approaches infinity. 

2. The Maximum Likelihood estimate is asymptotically Gaussian. 

A 

3. The Maximum Likelihood is invariant, i.e., if 0 is the Maximum 
Likelihood estimate of the parameter vector 0, and if u(9) is a 
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function of 6 with a single valued inverse, then u(0) is the 
Maximum Likelihood estimate of u(9) . 

4. The Maximum Likelihood estimate has a variance that approaches 
the Cramer-Rao lower bound asymptotically, i.e., it is asymptotically 
efficient. 

It can be easily shown that((17) and (19)), from a Bayesian point 
of view, the Maximum Likelihood parameter estimates for our model 
are unbiased. 

The next question of interest would be a comparison of the 
covariance of the parameter estimates obtained from different 
estimation techniques. It can be seen that an efficient estimator 
(i.e., that estimation scheme which gives the parameter covariance 
equal to the Cramer-Rao lower bound) can onl/ be the Maximum Likelihood 
estimator. An efficient estimator also has to satisfy 



where 0(Z) is an unbiased estimate of 0 (i.e., E(0(Z)) = 0), J is 

the likelihood function and k(0) is any function of 0. Hence, if 
the likelihood function J does not satisfy equation 34, then 
nothing is known about the covariance of the parameter estimates 
obtained by the Maximum Likelihood method. In the above condition, 
unbiased estimators which give lower covariance than the Maximum 
Likelihood estimator may exist, though there does not exist any 
general rule for finding them. 
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2.8 OUTPUT ERROR METHOD USING QUASILINEARIZATION 

We use an iterative method beginning with an initial parameter 

A 

estimate 0 * 0^. The problem is to find a zero of the gradient 
of the cost function 33, 3J/30 - 0. Consider a two-term Taylor 
series expansion of 3J/30 about the kth iteration value of 0 

(3J/30) k+1 Z (3J/ 30) k + C» 2 J/30 2 ) k A 0 k+1 (35) 

where 

A 

A ®k+l * e k+l * ®k 

2 2 

(3 J/ 30 is the second gradient of the cost function with respect 
to 0 at the kth iteration. If equation 35 is a sufficiently 
close approximation, the change in 0 for the (k+1) th iteration 
to make (3J/30) k+ j approximately zero is 

A ©k+1 * - C0 2 J/ae 2 ) k ]‘ 1 ( 3J/30 ) k ( 3 6)"^ 

Using now for v(j) the two term expansion 

V(j) k Z v(j) k-1 + (v(j)) k A 0 k (37) 


one obtains for the first and second gradients of the cost function 


M 


( 3 J/3 0 ) k = 2 


(3 2 J/30 2 ) k s 


j = 1 R ' 1 Cv(j)l k 

2 a# R - 1 


(38) 


(39) 
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£ 

We thus need solutions for v(j) k and -gg- vCj)^ . 
we first solve the system and measurement equations 


For this purpose 


1 

X * f ( X * u,t) 


(AO) 


y « h(x,u,t) 

for each iteration whereby the initial conditions are either obtained 
from the measurements or are included in the unknown parameters 9. 
The innovation is now obtained from equation 23. Next we solve the 

"sensitivity equations" for each iteration 

i 

ix/80£ * af/J)0£ ♦ 3f/3x I * »x/a0£ 


3y/&0£ * ah ^ ax |x«x ax / 30 i 


(41) 


A 

The initial conditions of Ux/39^ are zero except when x(0) is 

identified as part of the parameters 9. In this case the initial 

partials have the value one. With equation 23 we can now compute 

the first and second gradient of the cost function, equations 38 

39, and then obtain the change in parameters for the next iteration 

from equation 36. This involves the inversion of the sensitivity 

matrix M (equation39) , whereby M _1 is the Cramer-Rao lower bound 

for the covariance of the parameters. 

The method is easily extended to the case with a priori informa 

tion on th*» parameters, equation 25. The sensitivity matrix 39 is 

then augmented by the term 2P _1 , and the gradient 38 is augmented 

6 

by the term 2P" 1 (9-9 ), see (16). 

9 0 k-1 
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2.9 PARAMETER ESTIMATION BY FILTERING 

The parameter estimation methods discussed so far can be denoted 
as "global" methods. The performance criterion includes the test 
data for the entire duration of the transient. Filtering is an 
important tool in state and parameter estimation. It can be used 
either in conjunction with global estimates, or it can be used as 
a direct approach to state and parameter estimation. An example 
of the fir^t type of filter applications is the prefiltering of 
test data before using them in a least squares regression estimate, 
see for example (3). The Graham digital filter can remove high 
frequency noise. A Kalman filter can be used to estimate state 
variables and their rates not directly measured. It also removes 
the noise in the measurements. The role of the Kalman filter in 
Maximum Likelihood estimation has been shown in equations 27 to 31, 
where it is used to establish the innovation sequence. 

In addition to applications in global estimation methods, 
filters can also be used as substitutes for global methods. The 
advantage of such direct filter methods is a reduction in computer 
effort particularly in cases with a large number of parameters. The 
disadvantage is that unlike the inverted information matrix of the 
Maximum Likelihood method that provides a lower bound on the parameter 
covariances, no physically meaningful parameter covariances are 
obtained with the direct filter methods. The covariance propagation 
equations require initial values that are usually impossible to 
obtain in any rational way. Though improvements of the filter solution 
(forward time integration) can be achieved by smoothing (backward 
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time integration), the final parameter covariances remain arbitrary, 
since they evolve from arbitrary initial covariance estimates. 

Assuming that all state variables and their rates have been 
either measured or are otherwise known from manipulating the measure- 
ment data, the unknown parameters, if they occur in linear form in 
the state equation, can be found by application of a linear filter, 
see for example (5) . The classical regression method is a special 
case of this direct filtering method, namely for infinite initial 
parameter covariances. In classical regression one obtains a single 
value of the error covariance matrix. The direct filter application 
allows the use of a finite initial error covariance matrix and it 
gives the evolution of this matrix as a function of time. One thus 
obtains an indication when to stop processing the test data after 
their information contents has been exhausted. As mentioned before, 
the absolute values of the error covariances are meaningless, since 
one usually does not have a rational way of establishing initial 
values for the parameter covariances. 

A method that appears to be economic of computer time for 
large numbers of unknown parameters was used in (3) for application 
to helicopters. The method consists of a simultaneous identification 
of states and parameters with the help of a non-linear filter. In 
other words, the unknown parameters are treated as additional state 
variables. Since there occur products of state variables and para- 
meters, the system equation is a non-linear one. The so called extended 
Kalman filter appears to be particularly useful for this purpose. 
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Either non-linear filtering alone or linear filtering in combination 
with smoothing is performed. The absolute values of the parameter 
covariances are again of no physical significance since they depend 
on the arbitrary initial values. In the following, a brief discussion 
of the direct use of filters in parameter estimation is given. 

2.10 LINEAR FILTER METHOD OF PARAMETER ESTIMATION 
2.10.1 Linear Sequential and Global Estimators 

In (5) the parameter identification is performed from a 
"system equation" 

6=0 (42) 

and a "measurement equation" 

G = h(x,0)+v (43) 

equation 43 is actually the system equation arranged in a form where 
the left hand side contains all terms that are free of the unknown 
parameters 0. If the system equation is linear in the state 
variables x and in the unknown parameters 0, h(x,0) is a linear 
function of the parameters. The noise vector v refers only to 
the terms on the left hand side of equation 43. The state variables 
that are multiplied by the unknown parameters in h(x,0) must be 
noise free. To obtain the parameter 0, both c and x must be 
known. If only part of the variables in ? and x have been measured, 
Kalman filtering is required in order to reconstitute the missing 


terms. 
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Optimal parameter estimates 0 can be obtained, under the 
assumption that v is zero mean Gaussian white noise, by minimizing 
the cost function 

T 

j*(i/2){(0(o)-0(o)) T p: 1 (o)(e(o)-e(o))+ m-hU.erfiT 1 

9 («> 

(c-h( x,0))dt) 

where the a priori estimates 0(0), Pq(0) are assumed to be given 
together with the noise covariance matrix R. The differential 
equations associated with this optimal problem are (see for example 
( 21 ) ) 

9 = P 0 (3h/38) T R" 1 (C - h(x ,0) ) (45) 

Pq = - Pq( 3h/30 ) T R“ 1 (3h/30)P e (46) 

These equations can be integrated with the aid of the initial a 
priori estimate for the parameters 9(0) and their covariance 
matrix Pq(0), which results in the optimal parameter estimate at 
each time t given the preceding measurements. Since the initial 
parameter covariance is usually not known and the assumed values 
are rather arbitrary, the matrix Pq from the integration of 
equation 46 is not a useful measure of the actual parameter covari- 
ance. However, once Pq has approached zero, the effect of any 
further measurements on the estimate 0 also approaches zero as is 
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evident from equation 45. Pg, therefore, is valuable in judging 
for what length of time the data should be processed. Equations 
45 and 46 represent the "linear estimator" used in (22) and (5) . 

Instead of the sequential estimation by integrating equations 
45 and 46 with some initial estimates 9(0) and Pg(0), one can 
also obtain a "global" estimate directly from equation 44. If one 

A 

assumes that one and the same parameter estimate 9 is valid 
throughout the time range from 0 to T, one obtains by setting 
3J/39 = 0 


* 

= [P~^0)+ f Oh/30) T R" 1 (8h/89)dt]" 1 [P“lo)0(O) 


Oh/ae )^” 1 (c - h(x,G)) dt] 


(See for example the appendix of (5) . ) A convenient assumption is 

P“*(0) = 0, which means an infinite initial parameter covariance 
9 

matrix. Then the above estimate, equation 47, reduces to the 
equation error estimate, equation 9. The initial estimate 9(0) 
is then not required and the evaluation of equation 47 is reduced to 
the determination of fixed boundary integrals, a matrix inversion 
and a matrix multiplication. The parameter covariance matrix at 
the time T is given by the first factor of equation 47 : 


p q (t) = [Plh 


Oh/3 0 ) T R* 1 ( 3h/30 )dt 3” 1 


(48) 
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which follows from the integration of equation 46, (see for example 
the appendix of (5)). Pg(T) from equation 48 can again be used to 
judge whether or not all the significant information contents has 
been extracted from the data. 

2.10.2 Iterative Equation Error Estimation with Updated Kalman Filter 
When using the parameter estimation methods of the preceding 
section, it is necessary to first determine, from the noisy deflection 
measurements, estimates for the deflections, for their rates, and for 
the accelerations. In (22) this was done by passing the noisy 
deflection data through a digital filter that takes out the noise 
above a certain frequency without distorting the signal in the low 
frequency range. The filtered deflections were then either 
differentiated twice, or a Kalman filter was applied in order to 
obtain the derivatives. Later studies in (22) showed large errors in 
the parameters for too low cut-off frequency of the digital filter. 

It was then decided in (22) to omit the digital filter and instead 
use the Kalman filter in an iterative way. In typical examples, 
it was found in (22) that the second iteration was as accurate as 
the result with the combined digital and Kalman filter. A block 
diagram of the method is shown in Figure 4. The iteration loop is 
indicated by the double lines. The system Kalman filter 
gives optimal state estimates from incomplete and noisy input 
and output measurements. The filter needs estimates of the system 
parameters that are updated after each iteration. 
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,rror Estimation with Updated Kalman Filter. 
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Only simulated noisy blade flapping measurements were used in 
the Kalman filter. The filter provided the deflection rates and 
accelerations needed for the "global" parameter estimate, but not 
the deflections themselves. In other words, the parameter estimate 
was performed with the simulated noisy deflection measurements and 
with the rates and accelerations from the Kalman filter. In the 
first iteration, a Kalman filter with estimated parameter values 
was used (typically 20% error). After updated parameter values 
had been obtained, a second pass with an updated Kalman filter was 
performed, etc. The deflection data remained the same for each 
iteration, but the rates of deflection and the accelerations were 
updated. This method worked well for the single blade identification. 
2.11 BAYESIAN ESTIMATION AS A FILTERING PROBLEM 

If we extend the quasi-Bayesian or maximum a posteriori 
probability (MAP) criterion 16 to include both the parameters 0 
and the states x(t), we have 


max p(x(t) ,0/Z) 'v max p(Z/x(t) ,6)p(x(t) ,0) (49) 

x(t),0 x(t),0 


Assuming now the non-linear system and measurement equations 1 
and 2, and assuming further that states and parameters have Gaussian 
distributions of the form of equation 18, the criterion 49 becomes 
(see (3) and (23) ) one of minimizing the quadratic function 


J 


( 1 / 2 ) 



0 ( 0 )- 0 ( 0 ) 



♦ 



+ 


T 



)-h(x,u,t ) 



l|w(t)|j 


2 

rQ -1 r 



( 50 ) 
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subject to the constraint equation 1. If the system and measurement 
equations 1 and 2 are linearized about the current estimates x and 0 
the recursive solution of the minimization problem 50 results in the 
extended Kalman filter equations given for the continuous case by 
(see (3) ) 

x * (3f /3x )x + (3f/3u)u + POh/ax) T R" 1 (y-(3h/3x)x) 

0 = P 0x Oh/3x) T R~ 1 (y - (3h/3x)x) 

p * (3f/3x)P + P(3f/3x) T - P(3h/3x) T R“ 1 Oh/3x)P + Q (51) 

+ (3f/30)P 0x ♦ COf/30)P ex ] T 
Pqx - P 0x (3f/3x) T + P* ( 3f/30) T - P^Oh/Sx^R-J-Oh/SxjP 

with P» = - P Qx (3h/3x) T R _1 (3h/3x)P 0x 

Even if the original system is linear, the augmented system is 
non-linear and hence the filtering problem must be solved by a 
non-linear filtering technique. In (3) the raw data are pre- 
processed by a digital filter and by a Kalman filter that does not 
use the unknown parameters but merely makes use of the transformation 
equations from a space- fixed to a body- fixed reference system 
(Euler equations) . 

Lebacqz in (24) applies basically the same method except for 
a discrete instead of the continuous filter formulation. He further 
uses a one stage filtering- smoothing algorithm which has the advantages 
of reducing the bias due to non-linearities and of making the algorithm 
less sensitive to initial conditions. Mehra, in (1), is critical 
of using an extended Kalman filter for the augmented state including 
the unknown parameters. His arguments are that the uncertainties 
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in the states are usually much smaller than the uncertainties in the 

parameters. Therefore the assumption of local linearization about the 
latest estimate which are acceptable for state estimation with an 

extended Kalman filter are generally less valid for parameter 
estimation. Moreover, the filter for the augmented state assumes 
knowledge of the a priori parameter covariances which are unknown. 

As mentioned before, the arbitrary a priori parameter covariance 
used as initial conditions for a filter that includes parameters 
as state variables gives unreliable confidence limits on the para- 
meter estimates. An added difficulty of applying a filter to the 
augmented state is that poor a priori estimates of the parameters 
make the convergence rate slow or may even cause divergence of the 
filter solution. Though improvements can be applied to the extended 
Kalman filter like local smoothing and local iteration and smoothing, 
the basic shortcomings of this method appear to have been correctly 
described in (1). Unfortunately, the application of the complete 
algorithm of Maximum Likelihood identification given in (1) is for 
a large system much more demanding of computer size and time than the 
filter solution with the augmented state. While aircraft parameter 
identification with the complete Maximum Likelihood algorithm of 
(1) has not as yet been accomplished, the method of filtering the 
augmented state has been applied to several aircraft parameter 
identification cases, for example in (3) and (24) . 

2.12 IDENTIFIABILITY PROBLEMS 

Identifiability problems can occur no matter what identification 
algorithm is used. They are related to the initial 3 steps involved 
in system identification as listed at the beginning of this chapter: 



35 


the selection of a suitable input, the selection of the instrumen- 
tation, and the selection of the mathematical model. A few comments 
are added here to point out some difficulties that have been 
encountered due to these three initial steps. 

If the input does not adequately excite some of the system 
modes, the associated parameters cannot be adequately identified. 
Sometimes it is practical to combine the responses to various types 
of inputs into a single identification run, see (3) . While each 
of the single inputs excites only a limited number of modes, the 
combination of inputs provides an adequate excitation of all modes 
required for the estimation of the parameters. Efforts have also 
been made to design inputs on the basis of certain optimization 
criteria. More details on this problem are given in (25). 

If there are large unaccounted for instrumentation errors, 
non-physical parameter values may result. In (26), instrumentation 
lags and control measurement errors were found to be most significant. 
Static measurement errors and instrumentation lags can be a much 
greater source of parameter inaccuracies than white noise. A 
detailed analysis of the relationship between static and dynamic 
measurement errors in states and control inputs and the accuracy 
of the parameter estimates is required. 

If the selected mathematical model for the system is inadequate 
the parameters are forced to account for some unmodeled effects. 

The estimated parameters may, therefore, be quite different from 
those determined by aerodynamic theory or wind tunnel tests would 
indicate. A good example is given in (27) where a six degree of 
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freedom mathematical model for a helicopter gave unrealistic 
derivatives, since it had to account for effects of some neglected 
modes. A unique six degree of freedom linear model for the heli- 
copter flight dynamics does not actually exist. When a nine degree 
of freedom mathematical model is used, these difficulties disappear. 
Modeling errors are also a major cause for the lack of convergence of 
iteration procedures or of parameter identification by filtering 
methods. The best remedy against difficulties from modeling errors 
is the adoption of a more suitable mathematical model. Some other 
measures to improve the convergence of iteration procedures or of 
filtering methods will be briefly discussed. In the cases where a 
priori values of parameters, for example from theory or from wind 
tunnel tests, are available, one can use an a priori weighting 
matrix that expresses the confidence in these values and prevents 
the algorithm from deviating too much from the a priori values. 
Sometimes there exist some relationships between the parameters. 

These should then be used as constraint. 0 in the optimization problem 
to avoid non-physical parameter estimates. If parameter dependencies 
exist, difficulties are encountered in inverting the information 
matrix. An exact dependency between parameters should result in a 
zero eigenvalue of the information matrix. A rank deficient solution 
makes use of the fact that in case of near parameter dependencies 
there is a large spread between a set of small eigenvalues and another 
set of much larger eigenvalues of the information matrix. 

In filter solutions, divergence because of modeling errors can 
occur when the covariance matrix becomes prematurely too small, thus 
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prevonting further test data to be of influence. There are several 
ways to prevent premature small covariances. One can provide 
fictitious noise input to the system or one can directly increase the 
parameter covariance in each time step according to some rule. One 
can also overweigh the most recent data thus causing the filter to 
reduce its memory of the data of the more distant past. This 
indirectly increases the parameter covariance matrix. Since too 
short data length and too large errors in the initial parameter 
estimates may also result in non-physical parameter values or in 
divergence of the identification algorithm, longer transients and 
better a priori parameter estimates can lead to the avoidance of 
these difficulties. 

2.13 VALIDATION OF ESTIMATES 

Once a set of parameter estimates has been obtained the question 
arises: what confidence can be associated with this set? As 

mentioned before, the parameter covariance matrix obtained by 
filtering the augmented state is not a good measure of this confi- 
dence. The inverted information matrix obtained with the Maximum 
Likelihood method represents the Cramer-Rao lower bound for the 
parameter covariances and is a better measure of this confidence. 

Using the parameter estimates to predict the transients from 

which the estimates have been obtained, and computing the rms 

error with respect to the measured transients, gives another confidence 

measure. However, if the system is inadequately modeled, one may 

obtain a small rms error despite the fact that the parameter 

values are wrong in comparison to theoretical or wind tunnel results, (27). 
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A better way of validation is to compare the prediction with the 
results of test data not used in the identification process. In 
fact, it is good practice not to use all of the available test 
data for the parameter identification but to reserve some of the 
runs for such a comparison. Sometimes it is desirable to perform 
the parameter identification not just with one mathematical model 
but with a variety of models. In the case described in (27), 
a mathematical model with more parameters gave a much better identi- 
fication result than a model with fewer parameters, better in the 
sense of an improved correlation with theoretical}- and wind tunnel 
generated parameters. However, there are also cases where mathematical 
models with a larger number of parameters gave worse identification 
results than a model with fewer parameters, see (28). Adequate 
parameter estimation from transients requires careful attention to 
the many contributing factors in the input, instrumentation , 
mathematical modeling, and the estimation algorithm, and the 
validation of this process i only be considered complete after 
the rms errors of the prediction with the estimated parameters 
as compared to test data have been found acceptably small for all 
types cf possible transient excitations of the system. 

2.14 APPLICATIONS TO LIFTING ROTORS 

Lifting rotor characteristics are not well approximated by the 
usual set of aerodynamic derivatives. One reason is blade modes 
that must be considered particularly in rapid transients. Another 
reason is the dynamic rotor wake that is produced by the time 
varying rotor thrust . nd rotor pitching and rolling moments and that 
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has a feedback effect on the rotor forces and moments. The omission 
of the blade modes, as shown in (27), results in non-unique and 
non-physical rotorcraft derivatives. The identification is better 
if separate rotor degrees of freedom are introduced even in the 
crud» form of a first order lag as was done in (4) . 

A variety of identification methods has been used with respect 
to lifting rotors. After preprocessing the test data with a 
digital filter followed by a Kalman filter that does not contain the 
aerodynamic derivatives (transformation or Euler equations), least 
squares identification is applied to rotorcraft transient flight 
test data in (3) and (27). Each identification run is made with 
several transients simultaneously. The least squares results are 
then used as start-up values for an extended Kalman filter for 
the augmented state. It is not obvious that the extended Kalman 
filter actually improves on the least squares results, though filter 
convergence is achieved. In (4) the output error method with quasi- 
linearization is applied without preprocessing the flight test data. 
Hie flight data of both (3) and (4) were obtained in calm air. The 
equation error method in its filter form was applied in (5) to 
simulated noisy blade flapping and torsion measurements at high rotor 
advance ratio. The simulated data were preprocessed by a Graham 
digital filter, but not by a Kalman filter. (5) assumed that all 
states and their derivatives had been measured. In contrast (22) 
assumed that only flapping deflections are measured but not flapping 
rates or flapping accelerations. For the dynamic wind tunnel tests 
simulated in (22) there is no way of applying a Kalman filter that 
does not contain the unknown parameters. However, it was found in 
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(22) that for the cases studied, a Kalman filter with considerable 
errors in the unknown parameters was useful in obtaining the non- 
measured flapping rates and accelerations. The parameter identifi- 
cation was then performed by the equation error method in its filter 
form. 

In (29) the same method (except for using global estimates) is 
used in an iterative form. In addition, the c itput error method with 
quasi linearization is applied to the same and to more complex rotor 
identification problems. 
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3. MATHEMATICAL MODELS OF TOE FLAPPING RESPONSE 
OF A HINGELESS ROTOR 


3.1 SINGLE BLADE MODEL 

Using the simplest analytical model of a lifting rotor, a 
straight blade flapping about the rotor center, one has in a rotating 
frame of reference for the flapping "ngle 3 the following equation (15) . 

B + (y/2)C(t)B + o>,2b + (y/2)K(t)3 = (y/2) (ne e + ^ A) (52) 

One rotor revolution corresponds to t = 2ir. For neglected reversed 
flow effects, zero root cut-out and with tip loss factor B, the 
functions C(t), K(t) , mg(t), m A (t) in terms of rotor advance 
ratio p are (15) : 

C(t) = (1/4)B 4 ♦ (1/3 )B 3 p sin(t) (53) 

K(t) *<1/3 )b 3 p cos ( t ) ♦ (1/4)B 2 p 2 sin(2t) (54) 

m 0 (t) = ( 1/4 ) ( B 4 + B 2 p2) + ( 2/3 )B 3 p sin(t) 

(55) 

- (l/4)B*p 2 cos ( 2t ) 

m A (t) * ( 1/3 )B 3 + (1/2 )B 2 u sin(t) (56) 

In the numerical analysis, we use B = 0.97. A simple improve- 
ment of this analytical model that takes into account blade bending 
flexibility is possible (30). In transient conditions, the inflow 
A includes the dynamic rotor wake in a complicated form. 

As a first approximation of dynamic rotor wake effects one can 
use in equation 52, instead of the actual blade Lock number, an 




equivalent smaller value of y as suggested in (8) and (13) . 

Such an approximation can be expected to be satisfactory if the 
transient is relatively slow. For transients with high frequency 
contents, this approximation is invalid (11). In (11) it is seen 
that y*/y can be expressed by 

Y*/Y = 1 - 1/(1 + 8v/oa + 16Kji(»/oa) (57) 

The parameter v is defined in equation 70. 

The above formulation is based on single harmonic balance 
of the rotor root moment equation. The y* formulation reduces 
to the one obtained by momentum theory (9) (equation 58), when the 
phase variation is neglected and v = y (i.e., when X = v = 0) . 

Y*/y = 1/(1 ♦ ao/8y) (58) 

Due to rotor induced cross flow in a wind tunnel, the inflow 
parameter X will usually not be well known. In addition, the 
aerodynamic pitch angle 0 O is also not well known due to airfoil 
inaccuracies and pitch setting errors. For the wind tunnel tests 
considered here, we assume X = 0 and use the equivalent Lock 
number y* as an unknown parameter to be determined from the blade 
flapping measurements. In addition we have a transient blade pitch 
input 8b assumed to be known. The problem then is to determine 
from blade flapping transients caused by blade pitch inputs, the 
equivalent Lock number y* and the equivalent collective pitch 
setting 0 O . 

In order to obtain a more realistic description of the rotor 
dynamic inflow, it is necessary to formulate the rotor theory in 
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multiblade coordinates, as is done in the following section. The 
identification of y* is also possible with multiblade coordinates 
and has, for experimental data, the advantage that the measurements 
of all the blades are used, and hence, an averaging effect results. 

In fact, after it was found out that the blade pitch angles differed 
for individual blades by fractions of a degree, only multiblade 
measurements were used for the y* identification. 

3.2 MULTIBLADE FLAPPING EQUATIONS WITH DYNAMIC ROTOR WAKE 

It has been noted (10) , that for most purposes it would appear 
adequate to consider the first rotating mode elastic bending effects. 
The moment balance of all the flapping forces on a rotor blade 
about vhe hub is given by (see (10) ): 


(l/y)$ + (l/2)C(t)$ + (l/y)u>*8 + (l/2)K(t)8 = (1/2) (Am x + OBq) (59) 
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t 

x n dx 


(60) 


Here the actual first rotating mode n(x) is replaced by the closed 
form expression (see (10) ) : 

n » x ♦ icTsinhO. 93x)/2 sinh 3.93 ♦ sin(3.93x)/2 sin 3.93.1 

( 61 ) 

c. = 0 will correspond to a rigid blade mode. 

Since the dynamic rotor inflow that couples the motions ot tbe 
various blades is included, a multiblade representation is necessary. 

The relation between single blade and multiblade variables for a 


4-bladed rotor is: 
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Blade flapping angle: 



h 

®n 


= (1/4) 


6 d 


1 1 1 1 " 


•1 

2 cos t -2 sin t -2 cos t 2 sin t 


02 

2 sin t 2 cos t - 2 sin t -2 cos t 


03 

3-1 1-1 


»4_ 


Blade pitch angle: ®k = ®o ~ ®I s i n 'he. + ®II cos *he 


( 62 ) 


(63) 


Induced flow: = v 0 + Vj(r/R)cos tk + Vjj(r/R)sin ^ 

The variable represents differential coning for the 4 -b laded 
rotor, whereby one pair of opposing blades cones up, the other pair 
cones down. Though a linear distribution of the induced flow over 
the radius is defined in equation 63, this assumption is not required 
for the parameter identification process. Different inflow distri- 
butions merely produce different values in the identified parameters 
but do not change the form of the equations. 

3.2.1 Flapping Equations without Reverse Flow 

For low advance ratios the region of reverse flow is small. 

Since this region is concentrated near the hub of the rotor, the 
reverse flow affects the moment balance, and hence the flapping 
response, very slightly. Hence, for low advance ratios (generally 
acceptable for .4) the flapping equations are greatly simplified 
by neglecting the reverse flow effects without an appreciable error 
in the flapping response. 

Substituting the transformation equations 62 and 63 in the 
flapping equation 59, the multiblade representation of the flapping 
response is obtained. The limits of the integrals in equation 60 are 
from zero to B to take into account the tip loss factor. 



The closed form first mode expression 61 is substituted in equations 
60 which, in turn are introduced into the multiblade flapping 
equations. After some algebraic manipulation, the flapping equations 
are obtained as follows: 

6 0 + <B 4 Y/8)0 o + t*\ $ o + (B 3 Y|i/12)(3 n + X n ) 

+ (B 3 y/6)A o - (B 2 yu 2 /8) sin 2t 0 d - (y/2)k (-.028 B q 
+ .014 p 3 n - .457|i + .535p z B d > 

* (B 4 Y /8 + B 2 yv 2 /8)0 o - (B 3 yp/6)0 1 (64) 

B r + (bVsJBj + (o»2 - DBj + (B 4 y /8) (8 ix + 

+ ^ 6jj + (B 3 yw/6)B q + (B 2 yp 2 /16)6jj + (B 2 yw 2 /16) (sin 4t Bj - cos 4t 8^^.) 

"I 

- (B Yu/6)(sin 2t B d + cos 2t Bj) 

+ (y/2) k(.028(B x + B IX ) + (.886)p B Q 
+ (.268)p 2 sin 4t Bj + (.268)p 2 (1 - cos 4t)B IX 
+ (,028)u sin 2t B d - (.886)p cos 2t B d ) 

* (B 4 y/8 + B 2 yu 2 /16)8 xx + (B 2 yu 2 /16) (sin 4t 0^ - cos 4t 0^^) (65) 
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&H + (B 4 Y/8)B n + <«J - DBjj + (B 4 y/8)(-B i + A n ) 

-2 + (B 3 yu/6)B q + (bV/^Bj + (B 2 yu/4)A o 

- (B 2 YW^/16) (cos 4t $j+ sin 4t 0j.j) - (B\y/6)(sin 2t B d - cos 2t Bj) 
+ (Y/2) k((-.028)w S q + (.028)^ - Bj) 

+ (.268)y 2 (1 - cos 4t)B^ - (.268)y 2 sin 4t Bjj 

- (.028)y cos 2t 8 d - (.886)y sin 2t 0 d ) 

- (B 3 Yy/3)6 o - (B 4 y/ 8 + 3 B 2 Y y 2 /16)e i 

l 

- (B 2 Yy^/16)(cos 4t 0j + sin 4t 8jj) 

B d + (B 2 Y/8)§ d + p d _ (B 3 Y y/12)sin 2t (Bj + 2 B n + Aj) 

+ (B 3 Yy/12)cos 2t (Bjj - 2 Bj + A n ) 

- (B 2 Y y 2 /8)sin 2t B q - (y/ 2)< (-(.014)y sin 2t (Bj + Bjj) 

+ (.014)y cos 2t (Bjj - Bj) + (.535)y 2 sin 2t B q 

+ (.443)y cos 2t Bj + (.443)y sin 2t 8 IJ; - (.028)8 d ) 

= (B 2 Y y 2 /8)cos 2t e Q - (B 3 Y y/6)(cos 2t 6 X + sin 2t e n ) 


(67) 
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3.2.2 Inflow Model 

We adopt here the rotor inflow perturbation model of (9) and 
(11). The inflow model is based on the relation between the aero- 
dynamic rotor thrust and moment coefficients and the perturbation 
inflow state variables. Equation 33 of (11), written in our 
notation, reads 
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( 68 ) 


Rotor thrust and moment coefficients Op, Cj^, Cg are from aerodynamic 
contributions only. Lg is the empirical L-matrix defined in (9) . 

The theoretical values of k M , kj^ and kj^, using potential flow 
around a solid disk are given in (9) as k^ = .849, kj^ = kj^ = .113. 
The components of the L-matrix as well as kj^, kj^ and kj^ will be 
identified from rotor transient tests. From momentum theory, one 
obtains according to (11): 

r 

2v 0 0 

0 -v/2 0 | (69) 

0 0 -v/2 


-1 J_ 

[L] = oa 


with 


v - P 2 + X(X -f v) 

V ‘ (p 2 + P)l/2 


( 70 ) 


where X and v are the trim values, about which the rotor inflow 
perturbations v Q , Vj , Vjj are taken. 
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The complete matrix can be defined as a matrix consisting 

of 9 parameters as follows: 


-1 


Ji- 

ao 


L 11 

L 12 

L 13 

l 21 

L 22 

L 23 

L 31 

l 32 

L 33 


PI) 


Substituting equation 71 in equation 68, we get, on expanding 

( L 11> v o + CLl2)Vj ♦ (Li 3 )v n = (1/KjpC r (72) 

Vi - Ud/K^) (L 21 )v c ♦ (L 22 )v x + (L 23 )v n = -(l/Kj^Cj, (73) 

v„ - n(l/K l2 ) (L 31 )v d + (L 32 )vj ♦ (L 33 )v n = -d/K l2 )C L (74) 

The parameters t q , t j and tji are now defined as follows: 

T o * ‘/V T J ‘ '^i!' s 1«I 2 
The thrust and moment coefficients C T> and are obtained 
as a function of the state variables, the details of which is given 
in 8.2. 

To match the perturbation inflow model (equation 68) , where the 

inflow variables v , v,, and v TT are considered as perturbations 

oi II 

about its trim values v" o , and Vjj, the flapping equations have 
to be perturbations about certain trim inputs. Since the flapping 
equations 64, 65, 66 and 67 are linear equations and since C^, 

C in equations 72 to 74 are linearly related to the state variables 
( see section 3.2 ), the state variables in equations 64 to 67 and 
in equations 72 to 74 can be considered as perturbation variables. 
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X , X , X now are identical with the induced downwash variables 

O I 11 


Hence, X^ = v„; X T = v T and Xjj = Vjj. 


V V I* V 1I* » "o _ 'o’ "I ” V I 


3.2.3 Flapping Equations with Reverse Flow 

With increasing advance ratio, the region of reverse flow becomes 
larger and its effects can no longer be neglected. The limits of the 
integrals in equation 60 are no longer simply from 0 to B, but 
are split up depending on whether the flow in the region is normal, 
mixed or reversed. 


Region 1 : J * 

Region ? : I 


ft 


Region 3: 
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i. 

f - [ 

J-M sinif/ Jq 

f 


-U siniji 


The three regions are clearly explained in (15) . To obtain 
a closed form analytical solution for the coefficients of all the 
states in the flapping equations (and in the thrust and moment 
coefficients), a fourier expansion is obtained for all the coefficients 
around the azimuth. This gives rise to different sets of coefficients 
for different advance ratios. The flapping perturbation equations 
as well as the thrust and moment coefficients are obtained in a manner 
similar to those obtained by neglecting reverse flow. The blade is 
assumed to be rigid (i.e. k is assumed to be zero). The coefficients 
are provided in 8.4. The flapping equations are listed in 8.3. 


3.3 EXCITATION OF PITCH STIRRING TRANSIENTS 


For wind tunnel experiments with pitch stirring transients the 
initial state of the rotor will be given by prescribing the advance 
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ratio, the collective pitch angle, the rotor angle of attack and 
the cyclic control setting that will be zero longitudinal cyclic and 
1.5° lateral cyclic. 

The lifting rotor wind tunnel model described in (31) allows 
excitation of progressing and regressing flapping modes at various 
frequencies. By a minor modification of this model, progressing 
or regressing transients can be excited. One can describe such inputs 
as pitch stirring transients. In a helicopter, this would amount 
to cyclic stick stirring, whereby the amplitude of the cyclic pitch 
would remain constant while .he frequency of the stirring motion 
changes. At the time t Q , pitch stirring is initiated. If we 
denote the angular pitch stirring speed as u>, positive in the 
direction of rotor rotation, and the pitch stirring angular acceler- 
ation as <L, assumed to be constant, we have 

u = w(t - t Q ) (75) 

For a progressing mode <■> is negative and for a regressing 
mode u is positive. In a rotating reference system the blade pitch 
angle is given by 

0 = 0 O ♦ 1,5 cos t«(t-t 0 ) ♦ t] 


w 


0 for t < t 0 

«(t - t 0 ) for t > t c 


In a multiblade representation the blade pitch angle of the kth 
blade is 


(76) 


9. = 9 - 9, sin +0 cos * 

k o I k II k 


( 77 ) 
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where 


•l 


0 for t 4 t Q 

1.5 sin «o(t-t 0 ) for t > t 0 


®II 


1.5 for t 4 t 0 

1.5 cos u)(t-t Q ) for t > t 0 


(71) 


(79) 


The meaning of these input equations is the following. At the 

time t * 0, a step lateral cyclic pitch input of 1.5 degrees is imposed. 

At time t * t , the response to this input is approximately 
o 

e 

stabilized. At this time the pitch stirring acceleration of u> is 
introduced which leads to a progressing flapping excitation. The 
identification starts at t = t with the pitch stirring transient. 

9j represents forward cyclic pitch, 9 jj represents left cyclic 
pitch. If perturbation equations are considered, the perturbation 
at time t Q is zero, 9j excitation stays the same but 0 II excitation 
is now defined as: 



0 


for t < t 


0 


(» 0 ) 


I 1.5 cos wCt-t^) - 1.5 for t > 

The wind tunnel experiments are conducted with a variety of pitch 
stirring accelerations. Generally, the computer experiments are 
conducted with a pitch stirring acceleration of 


a) = - . 1/ir 

which is in the progressing sense. 


( 91 ) 
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Since in the non-dimensional time units used here the time of 

one rotor revolution is 2ir, the angular pitch stirring velocity one 

rotor revolution after initiation of pitch stirring is .2, that is 

one fifth of the rotor angular speed. Figure 5 shows the ti^e history 

of the blade pitch for about two rotor revolutions (t Q - 0, t * 0 

♦ 

to 12) in a rotating frame of reference for w = -. 1/v, corresponding 
to equation 76. Figure 6 shows the time history of blade pitch in 
multiblade representation, that is 0j and 0jj vs. time t for the 
same acceleration corresponding to equations 78 and 80. Figures 5 
and 6 refer to the progressing mode. The physical reason why 
regressing modes are less suited for rotor wake identification is the 
fact that at a certain regressive excitation frequency the excitation 
is in resonance with the regressing flapping mode. At this condition 
no induced dynamic rotor wake exists since aerodynamic excitation 
and aerodynamic damping cancel each other. Since regressing mode 
transients include a frequency region with a weak dynamic rotor wake, 
the identification of the wake parameters is not expected to be as 
good as it is for progressing mode transients. 
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Figure 5. Time^History of the Blade Pitch in a Rotating Frame of 
Reference for w • — . 1 /ir in Equation 76. 
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Figure 6. Time History of the Blade Pitch o and t) tor u ■ -.1/it 
in Equations 78 and 80. 1 
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4. SIMULATION STUDIES 

4.1 SELECTION OF IDENTIFICATION METHOD BASED ON PRELIMINARY 

SIMULATION STUDIES 

The selection of the identification method used in this thesis is 
based on several simulation studies ( (5) , (22) and (29) ) . 

Simultaneous state and parameter identification in (32) and (33) 
was conducted using an extended Kalman filter. A major drawback of 
the method is that the filter can easily diverge unless good initial 
estimates are available. Particularly in (33) a considerable effort 
was applied to obtain such good initial estimates. The test data 
was first processed with a digital filter that took out high frequency 
noise without distorting the main signals. The data was then processed 
with a Kalman filter based on the Euler equations, which do not 
contain the unknown parameters. Thus measurement bias was removed 
and missing cha nels were reconstituted. Finally a least squares 
algorithm was applied to obtain estimates of the unknown parameters. 

The subsequent application of the extended Kalman filter led to 
modified parameter estimates, however it is i. t clear whether or not 
these modifications represent improvements. In any case the modi- 
fications were not large, and the initial estimates appeared to be 
satisfactory approximations . 

In trying to apply the experience from (33) to wind tunnel 
model transients a difficulty arises, in that there is no equivalent 
to the Euler equations for the ai..raft. Thus there is no way of 
using a Kalman filter which is free of the unknown parameters. 

Instead, if a Kalman filter is to be applied, estimates of the parameters 
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must be inserted* Another difficulty for our wind tunnel model 
tests is that only flapping deflection measurements are made, while 
the rates of deflection and the accelerations are not measured* Thus 
the Kalman filter with the estimated parameters is called upon to 
provide both rates and accelerations. 

This method as explained before has several disadvantages. 

A priori parameter covariances being generally unknown, give poor 
and most often, wrong parameter covariances as solution of the 
Riccati equation. Wrong estimates of the parameter values often 
cause the filter to diverge. 

Analysis is also made by replacing the least squares algorithm 
of (33), by linear sequential estimators and a simpler "global" 
estimator. This has the advantage that finite initial parameter 
covariances can be used, and that the time history of the parameter 
covariance provides a measure for the time beyond which no more useful 
information can be extracted from the test data. 

The linear sequential estimator as shown before (used :'n (22) ) 
requires the simultaneous integration of the filter and of the covariance 
differential equations. A simpler "global” estimate requires only the 
inversion of a system of linear equations for the unknown parameters 
and the evaluation of a number of integrals over the time period of 
the transient. Therefore, a number of comparisons were made between 
these two methods. 

For a single blade 2 parameter identification , both the linear 
sequential estimator and the "global" estimator provide quite 
accurate parameter estimates. For convenience, > and 6 £ y instead 



57 


of y and 0 Q were identified. Preliminary analysis and comparison 
between the Iterated Equation Error estimation with updated Kalman 
filter and the Maximum Likelihood method had given the following 
results: For single blade parameter identification from pitch 

stirring transients the Equation Error method applied in an iterative 
form using a Kalman filter with the latest parameter updates worked 
well and required the least computer CPU time. For multiblade 
parameter identification this method became impractical because of 
slow convergence and high computer CPU time. The Maximum Likelihood 
method worked well both for single blade and multiblade applications, 
though in case of single blade identification it requires somewhat 
more computer CPU time. The parameter covariances from the Maximum 
Likelihood method are clearly superior to and more meaningful than 
the covariances determined with the Equation Error method. The 
Maximum Likelihood method also gave good parameter identifications 
in the presence of both measurement and system noise, though most of 
the computer experiments were conducted with measurement noise only. 

The following Table 2 compares the results of the various 
methods on the single blade model (29) . The last 4 rows refer to 
the Maximum Likelihood method. 

The number of iterations indicated in the table is that for 
which convergence has been achieved. The Iterated Equation Error 
estimation with updated Kalman filter needs the lowest total computer 
effort, however, the accuracy of the estimate is worst for y. The 
Maximum L"elihood estimation, due to faster convergence, needs only 
moderately more computer effort and yields better accuracy. 



Table 2. Comparison of Various Identification Methods on the Single Blade Model. 
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During the last decade the Maximum Likelihood method of para- 
meter identification has been successfully applied to airplane and 
helicopter transient testing. This method does not require 
preprocessing of the test data and also does not need complete 
measurements of the deflections, of their rates and of the accelerations. 
The parameter covariance estimates obtained with this method are more 
meaningful than those obtained with the linear sequential estimator 
used in (22) . 

From the above study, one can conclude that the Maximum Likelihood 
method in its simplified form in which system noise is not modeled, 
is, for the applications studied, superior to the Equation Error and 
other existing methods and thus will represent the method of choice 
for the parameter identification from wind tunnel rotor model tests. 

4.2 SIMULATION STUDIES FOR FORWARD FLIGHT USING MAXIMUM LIKELIHOOD METHOD 
The Maximum Livelihood method for our particular case pertains to 
the system equation (zero system noise) 

x = f (x, u, 0) (82) 

0 is the vector of unknown parameters that may include initial values 
of the state variables, constant measurement bias, etc. The measure- 
ment equation is assumed to be linear and of the form 

y = H x + v (83) 

y is the vector of observed quantities, H is a matrix relating the 
state variables to the observations, v is the vector of random 
measurement errors, assumed to be zero mean white noise with given 
covariance matrix R 

E[v(t )v T (t)] = R 6 (t , t) (84) 

R is assumed to be constant with time. Though the preceding 
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equations do not show bias terms, bias errors could easily be included 
in the unknown parameter vector 0. 

A sample of measurements . . . y n is now made during the 

a* 

time of the transient and the parameter estimate 0 is selected such 
that the conditional probability of this sample of measurements given 
0 is maximized. 

A 

0 a max p(y. ... y /0) (85) 

9 1 n 

The following steps lead to the maximum of the likelihood function 
p(yj . . . y n /9) , though there is no assurance that the maximum is 
global. The method outlined here is called quasi linearization with 
the modified Newton-Raphson method. It assumes Gaussian distributions 
of the random variables. 

1. Select an initial parameter estimate 0 = 9 Q . 

2. Solve the system equation 14 with this parameter estimate 

x = f(i, u, 0) < 86 > 

The initial conditions can either be obtained from the 
measurements, or, where this is not feasible, they can be 

A 

included in the unknown parameter vector 0. 

3. Calculate for each measurement the "innovation term" 

V j = Y i ' H <87) 

4. Solve the "sensitivity equations" 

3 x/e k = 3f/30 k ♦ F(t) 3 x/ae k (88) 


A 

X 


where 


F(t ) * 3 F ( t ) / 3 
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The initial conditions of dx/SO^ are zero except when 
x(0) is identified as part of the parameter vector 0. 

In this case the initial partials have the value one. 

5. The likelihood function for zero system noise is 

J = log p(y x . . y /0) * - (1/2) E vT R" 1 v (89) 

5=1 3 3 . 

Determine now the gradient of this function with respect 
to 6 

N T l 

3J/30 = E vt R 3v./30 (90) 

j = l 3 3 

where 3v^/30 * - H 3 x^/30 (91) 

6. Compute the information or sensitivity matrix 

N 

M = 3 2 J/30 2 = E (3v./30) T R" 1 3v./3© (92) 

5=1 3 3 

The inverse M“* of the information matrix provides a lower 
bound for the cc/ariance of the updated parameter estimates. 

7. The updated parameter estimate is 

0 = 0 + A 0 (93) 

o 

where A0 * -H** 3, (3J/30)^ ( 94 ) 

8. Go now back to equation 86 with the updated parameter 
estimate and repeat the steps to equation 93 . Reiterate 
until convergence of the information matrix and of the parameter 


vector is obtained. 
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The Maximum Likelihood method, which was used quite successfully 
in the single blade forward flight analysis and in multiblade hovering 
analysis (29) , was extended to the multiblade forward flight analysis 
with the time delayed rotor inflow L-matrix model as defined before. 

To study the question of L-matrix identification, simulation 
studies were performed for a hypothetical rotor with the 
characteristics of the model rotor treated in (9) , so that the 
experimental values of the L-matrix determined in (9) could be 
used. The rotor was assumed to have a blade solidity ratio of .100 
and a blade flapping frequency of o>j =1.20. The experimental 
model rotor has a blade solidity ratio of .154 and is usually run 
with a rotor speed corresponding to a blade flapping frequency of 
Uj = 1.17. An advance ratio of .4 is the upper limit for a conventional 
helicopter. The dynamic rotor wake effects are rather small at 
this rotor advance ratio. The dynamic rotor wake parameter identi- 
fication should be easier for lower advance ratios where the rotor 
wake has substantially larger effects. In order to study the feasibility 
of L-matrix identification, the most unfavorable case of .4 advance 
ratio was selected. 

The inflow model chosen is given by equation 63. The theoretical 
values of K^, and Kj^, using potential flow around a solid disc 
are = .849 and Kj^ = Kj^ A K j =.113. Choosing the values of 
’ the parameters of the L P matrix from (9) at p - .4, the inflow model 
is obtained as 
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Here the number 7.5 represents the theoretical value of the ratio 
Kj^/Kf. We will assume this value as given and identify only K^, 
but not and Kj separately. The inflow model is now assumed 
to have six unknown parameters 0j to 0^, so that the above equation 
95 is written in the form 


X 


~ e i 

0 

0 

V 


c T " 


+ 

0 

02 

©3 

V I 

= ©6 

- 7 - 5C m 



0 

0 4 

©5 

_ v n 


_" 7 - 5cl _ 


The flapping equations are the perturbation equations derived 
without reverse flow (equations 64, 65, 66 and 67). The system 
equations have the given parameters ao = 2ir/ 10 , u = 0.4, tip loss 
factor B = .97, and the flapping frequency coj = 1.20. The Lock 
number will be assumed given as y = 5.0 for some identification runs 
and will be assumed an additional unknown in other runs. 

The inflow equations represent a feedback system, whereby most 
of the unknown parameters occur in the feedback loop. The feed- 
back signal, v, is unknown. The only measured quantities are 
0j, 9 jj, 3 0 , Bj, Bji, and 3^. The four other state variables are also 
unknown. The identification problem thus has seven unknown parameters 
ii ( is included, four unknown state variables and three unknown 
feedback variables. One of the unknown parameters (9^) is a time 
constant. 

The following studies were made using the above inflow model 
(equation 95) 

a) The angular acceleration of the pitch stirring shaft was 

considered to be <I> = -.1/u. The time of transient measure- 

ment is t = 0 to 12 time units, the sampling rate 
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At * 0.1 and the standard deviation of the measurement noise 
o v * .05. Hie average least squares fit o^(fit) is also 
noted > to determine the accuracy of the identified fit (Table 3). 

Prom the above identification am it is seen that 
though most of the parameters converge, the accuracy of 
the estimates vary. The diagonal terms Gj, 83 and 85 of 
the L matrix converge to more accurate parameter estimates 
than the off-diagonal terms 0 j and 0 ^. y A 0 y is very 
sensitive to the response and hence converges faster and 
accurately to the true value. 

(b) A priori knowledge of the parameter estimates is added to 

the likelihood function as a quadratic term involving the 

weighted difference between the estimated parameter values 

and the priori parameter values , i . e . , 

(0i - 0 O ) T A (Oi - e o ) (97) 

where A is the weighting matrix. The identification runs 

were under the same conditions as in (a) . Two pitch stirring 

accelerations of -.05/ir and -.1/ir are studied. The weighting 

matrix A was taken to be 2001. The Lock number, y, which can 

be determined quite accurately, is assumed to be known at y s 5.0. 

The main point of interest in Tables 4a and 4b is the 
improvement in accuracy of M~* values for the different 
parameters from & = -.05/* to w = -.1/ir. The values in the 
latter are about one-tenth the corresponding values for the 
M“* values obtained by using the slower acceleration. 

There does not seem to be any appreciable change in o^(fit) 
between the two cases in Tables 4a and 4b and Table 3. 



Table 3. Simulation Identification Study Results for the L-matrix Model 

(Equation 95) using the Maximum Likelihood Method with w = -.1/ir; 
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Table 4a. Simulation Identification Study Results for the L-matrix Model 
(Equation 95) using the Maximum Likelihood Method with a Priori 
Knowledge of the Parameters ip form of (97) added to the cost: 
t = 0-12; At ■ 0. l;o = 0.05;w = -.05/*. 
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Table 4b. Simulation Identification Study Results for the L-matrix Model 
(Equation 95) using the Maximum Likelihood Method with a Priori 
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A priori weighting does not seen to have helped in Obtaining 
better estimates of the parameter values. In cases of para- 
meters 0§ and 0g, tiie identified values are much worse. Hie 
accuracy of most of the identified parameters don't seem to 
be very good. 

(c) Another study of interest is the case where the induced 
flow is represented by the equivalent Lock number, y*. 

Hie simulated measurement data is obtained by using the full 
L-matrix induced flow model in determining the flapping 
response. This response is polluted with zero mean, 

Gaussian white noise of o v * .05 to obtain the measure- 
ment data (as before). Hie identification results are 
given in Table 5. Other conditions are the same as in 
(a) and (b). 

Hie equivalent Lock number, y*, was identified and 
converges rapidly. Hie M~* value obtained is low which 
shows the goodness of the identified value; but incomplete 
modeling (i.e., lack of an inflow model) gives rise to a 
o^(fit) value which is approximately twice as much as 
that obtained by using a complete inflow model in Table 4. 

(d) In order to identify the elements of the Lg _1 matrix 
directly and also to identify the three mass and inertia 
terms separately, the inflow model as given by equation 68 
is used rather than the simplified equation 96. 
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Table 5. Equivalent Lock Number y Identified from Data 

Generated from the Full L-matrix Induced Flow Model 
t * 0 - 12; At * 0.1; o v * .05; & s -.Oyv. 
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( 98 ) 


where: 0 6 A 1/K^ 

0 7 A 1/K 

l l 

0 a A !/ K t 

8 * h 

The identification program was run under the same 
conditions as in (a) except that length of the identified 
data length is t = 0 to 18 instead of t = 0 to 12. 

The increase in data length was justified on account of 
the poor accuracy of the identified parameters and the 
high values of the Cramer-Rao lower bounds for the 
parameters. This increase in data length was verified in 
a later chapter on "Optimal Data Utilization". The 
results are tabulated in Table 6. 

Once again, the off-diagonal terms of the Lg - * matrix, 
show poor identifiability. The (fit) value does not 
show noticeable change from the results in (a) . The improve- 
ment in the Cramer-Rao lower bounds with increase of the 
data length was expected. 

The simulated measured responses together with the 
identified responses are given in Figures 7, 8, 9 and 10. 

The perturbation identified inflow is plotted in Figures 
11, 12 and 13. 



71 


Detailed analysis of the simulation studies for pitch 
stirring acceleration of & * .05/w, -.1/ir and .2/v 
are given in (34). 



Table 6. Simulation Identification Study Results for the L-matrijc Model 

(Equation 98) using the Maximum Likelihood Method with w ■ >.l/v; 
t ■ 0 - 18; At ■ 0.1;e » 0.05. 
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4.3 EIGENVALUE ANALYSIS 

Hid non-uniform downwash is strongly coupled to the moment 
response of hingeless rotors (9) and (11) . This led to the study 
of the effect of the unsteady downwash on the rotor transient response 
and also rotor stability. Study of rotor transient response is made 
in the identification analysis of dynamic rotor response in the 
previous chapter. Stability studies are made by eigenvalue analysis 
of the dynamic flapping response of a rotor model. 

As mentioned before, it was seen in (14) that, in hovering, 
the damping of the regressing flapping mode is substantially reduced 
by dynamic inflow effects at low collective pitch angle. In the 
following, results of the forward flight eigenvalue analyses are 
presented using different analytical models. 

A detailed study was conducted using the forward flight model 
given by equations 64, 65, 66 and 67. The parameters of the inflow 
model given by equation 68 is obtained from Figure 4 of (9). Several 
important aspects of this analysis are: 

The eigenvalue analysis is first conducted at n = 0.4 using 
a complete flapping transient model (including the feedback downwash 
model) with and without the periodic terms. The comparison of the 
eigenvalues is given in Table 7 . Another interesting comparison 
with the above are the eigenvalues obtained by neglecting the down- 
wash in the constant model. For the constant system, the eigenvalues 
are obtained directly by taking the Laplace transform of the system 
equation and then solving for the roots of the characteristic equation. 
The periodic system equations on the other hand, have to be solved by 



81 


Floquet theory as outlined in (35). According to the Floquet theory, 
the imaginary parts of the eigenvalues are indeterminate and in their 
multiples of one can be added or subtracted. 

From Table 7 one can see that the error in representing the 
periodic system by a constant system (obtained by neglecting the periodic 
terms) is small, and therefore, the constant system is in our case, an 
accurate representation for the set of parameters used. The error 
in the eigenvalues in the above comparison are all within 2%. In 
contrast, the flapping model without the dynamic inflow model is 
significantly in error. The dynamic inflow reduces the damping of 
the regressing mode (fourth row of Table 7) by 40%. 

The sensitivity of the eigenvalues with variation in parameter 
values was also studied for the constant model. The flapping eigenvalues 
seem to be insensitive to changes in values of the L-matrix parameters 
in the downwash model. Variation in value of the Lock number y 
caused the real part of the eigenvalues to move closer to the imaginary 
axis with decrease in the value of y. 

The sensitivity of the eigenvalues to the parameters y (Lock 
number) and (blade natural frequency) is determined. The model 

chosen was the constant system without downwash. Change in a>i only 
changes the frequency component of the eigenvalues. The real part 
of the eigenvalue stays steady at A/2 (i.e., approximately y/16). Table 


8 shows the details. 
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Table 7. Comparison of the Eigenvalues Between the Three Different 
Forward Flight Models of the Flapping Response. 


Eigenvalues at y = 0.4 (y = 5.0; wj 2 = 1.4) 


Periodic System 

Constant System 
with Downwash 

Constant System 
without Downwash 

-0.25610. 137j 

-0. 25312. 134j 

-0.27412. 16 j 

-0.24510. 131 j 

-0.24611. 130 j 

-0.27611. 161 j 

-0. 27511. 162j 

-0. 27711. 168j 

-0.277+1. 168j 

-0.20010. 190 j 

-0.201+0. 194j 

-0.280i0.167j 

-0. 681+0. Oj 

-0. 682+0. Oj 


-1. 299+0. Oj 

-1. 316+0. Oj 


-1. 682+0. Oj 

-1. 665+0. Oj 











83 











84 


The variation in the eigenvalues from hover to an advance ratio 
of )i * 0.4 is given in Table 9. The constant system model with 
the complete downwash model is used. The parameter values used 
are those at p * 0.4, and given in equation 95. According to (9) 
the dynamic inflow parameters do not change much between p = 0 and 
P - .4. 

Surprisingly, there is negligible change in the eigenvalues 
with change in advance ratio. The variation may he more pronounced 
if the changes of '•'he inflow parameter values with change in advance 
ratio were taken into account. 

At higher advance ratio the effect of reverse flow and periodic 
terms becomes important. To study this effect a comparison between 
three cases is made as shown in Table 10 for p = 0.8. Eigenvalues 
are collared for the following three models: 

(1) A rotor model with the periodic terms, reverse flow 
effects and the complete inflow model 

(2) The above complete model neglecting the reverse flow 
and periodic terms 

(3) The model given in (1) with the downwash equations 
neglected. 

From Table 10 it is seen that, neglecting the reverse flow and 
periodic terms did not affect the flapping eigenvalues significantly. 
The eigenvalues corresponding to the downwash were changed greatly. 
At high advance ratio, the feedback due to the dynamic inflow 
becomes relatively unimportant which is clearly seen in Table 10. 
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Table 9. Sensitivity of the Eigenvalues to Variation in Advance 
Ratio in a Constant Forward Flight Mathematical Model 
with Downwash. 


Constant System with Downwash 


. 25912. 143j; -0.27711. 168 j; -C 

>. 25711. 157j; 

. 717+0j j - 1. 245+0 j; -1.673+0js 

1 -0. 19710. 197J. 

.25912. 142 j ; -0.27711. 168j; -t 

1.25611. 137 j; 

. 715+Oj ; -1.2S+0j; -1.672+0j; 

-0. 19810. 197j. 

.25812. 140j; -0.277±1.168j; -t 

L 25411. 136 j; 

. 708+0 j; -1. 264+0 j; -1.671+0jj 

; -0. 19810. 196J. 

.256±2.138j; -0.277±1.168j; -C 

1. 25111. 133 j; 

.697+0j; -1.286+0j; -1. 669+0 j ; 

; -0.19910. 196 j. 

.253±2. 134 j; - 0.27711. 168 j; -C 

1.24611. 130j; 

.682+0j; -1. 316+0 j; -1.665+0j; 

; -0. 20110. 194J. 
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Table 10. Comparison of the Eigenvalues Becveen the Three Different 
Poniard Flight Models at High Ad/ance Ratio (ji = 0.8; 

« 1.4; Y * 3.2). 


Model (1) 

Model (2) 

Model (3) 

-.18112. 143 j 

-.16712. 143j 

-.18712. 148j 

-.18711. 126 j 

-.174±1.147j 

-.18711.1513 

-.18711. 147 j 

-.177±1.166j 

-.18711.1513 

-.176t0.190j 

-.180±0.175j 

-.187i0.154j 

-. 580+0 j 

-.268+0j 


-2.45110. 151j 

-5.148+0j 



-17.937+0j 
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4.4 OPTIMAL DATA UTILIZATION FOR PARAMETER IDENTIFICATION PROBLEMS 

WITH APPLICATION TO LIFTING ROTORS 

In aircraft or wind tunnel transient testing the question comes 
up as to what kind of transient should be selected. If the transient 
is too short, the parameters will be identified with inadequate 
accuracy. If the transient is too long, an unnecessary amount of 
data must be processed. The question we pose here for the Maximum 
Likelihood method is — given a required accuracy of the parameter 
estimate, and given an input function, what is the minimal quantity 
of measured data necessary to achieve this accuracy? There are some 
recent studies where certain criteria were used to define an optimum 
input. We will first briefly discuss two of these optimal input 
proposals, and then proceed to develop the method of optimal data 
utilization for a given type of input. 

4.4.1 Two Proposals for Optimal Input Design 

General questions of input design are: 

(a) What type of input function should be used? 

(b) For what time period should the response data be processed 
to enable identification of the system parameters with 

a specified accuracy? Are certain time periods of the 
response particularly rich in information contents and 
should they, therefore, be preferably used? 

There usually are some constraints on the input design like amplitude 
constraints, smoothness constraints (step or impulse inputs are 
mathematical idealizations but often practically not realizable). 
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instrumentation constraints, and constraints imposed by the selected 
analytical model that usually filters out the higher frequency 
contents of the input. 

Analytical solutions of the problem of optimal input design 
require the minimization of a cost function. Stepner and Mehra 
(1) use the sensitivity of the system response to the unknown 
parameters as the performance criterion for optimal input design. 

The time of the transient is assumed to be fixed. Thus questions 
(b) are not involved. The measurement equation is 

y m (t) = y(x, 9, u, t) + v(t) ( 99 ) 


We write the Taylor expansion with respect to the parameter 9 about 
the a priori estimate 9 Q of 9 and neglect higher order terms : 


y m (t) = y(x * 0 °* u » t > + Ye - 


yU» 0 O . u, t)(e-e 0 ) 


( 100 ) 


+ v( t ) 


In the output error method (9 - 9 0 ) is determined by a least 
squares solution of equation 100 for a fixed time period (t 0 , t£) . 
For a high degree of accuracy in determining (9 - 9 Q ) the sensi- 
tivity function 3y/39 must be large. The scalar performance index 
selected in (1) is 

J = Trace (WM) (101) 

where 
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Due to the introduction of IT* in M, the performance criterion 
favors the measurements which are more accurate. The weighting 
matrix W is based on the relative importance of the parameter 
accuracies. 

If we assume linear system and measurement equations 



(103) 

(104) 

(105) 


the optimum input u can be determined as a two point boundary 


value problem whereby the Hamiltonian includes the term 


(1/2)Mq (u T u - E/t f ) . 


The scalar y G is the time invariant Lagrange factor to be evaluated 
from the Euler differential equations of the optimization problem. 

It should be noted that the " energy constraint" equation 105 has 
no physical significance but is a convenient device to obtain smooth 
input functions. Physically, the input will usually be limited 
by amplitude rather than by the quadratic criterion (equation 105) 
and quite different "optimal" inputs can then be expected. 

(36) attacks the problem of optimal input design in an entirely 


different way as a time-optimal control problem by minimizing 



( 106 ) 
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under the following constraints: 

System Equations 

x * f(x, u, 0, t) , X(t 0 ) = x Q 

Sensitivity equations 

3x/3x c = (3f/3x)(3x/3x 0 ) , 3x(t Q )/3x 0 = I 

3x/30 = (3f/3x)(3x/36 ) + 3f/30, 3x(t Q )/3e = 0 

Information matrix equations 

M* 1 = -M” 1 (3v/30) T r" 1 (3v/30)M" 1 
where v is the innovation: 

v s y m - y 

and where the information matrix M is given by equation 102. 

Finally, Chen assumes an amplitude constraint 

|u| 4 V (112) 

and he prescribes the trace of the information matrix for time t£ 

c ii (t f ) = a \ (ll3) - 

One can show that for linear input u(t) into the system equation 
and for an input matrix independent of any unknown parameter, the 
optimal input is of the "bang-bang" form between the amplitude 
constraints. The solution of this problem requires a computer 
search which was not performed (36). Rather, an arbitrary set 
of bang-bang inputs in the form of Walsh functions was shown to 


007 ) 

(108) 

(109) 

(110) 
t ( 111 ) 
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result in a specific case in lower values of M"*(tf) (given tf) 
than those obtained by using Mehra's "optimal input". This apparent 
contradiction can be explained by the differential equation 110 
governing M~*. For a particular value of M - * the rate of decrease 
of M"* with time is dependent on all elements of 

Oy/ae) T R“ 1 Oy/ae) 

while Mehra, in his criterion (equation 101) optimizes only the trace 
of WM. 

While the input amplitude constraint (equation 112) used by 
Chen is physically more significant than the quadratic constraint 
(equation 105) used by Mehra, the actual constraints are usually 
still more complex. In cases of airplanes or lifting rotors one 
usually wishes to limit the response to the linear sub-stall 
regime, since the analytical model to be identified is often a 
linear one. The stall boundary is, however, a complex function of the 
input and cannot be represented by an amplitude constraint for the 
input transient. This is particularly relevant for the lifting 
rotor, so that neither the Mehra nor the Chen input optimization 
criteria is useful for lifting rotor applications, quite apart from 
the excessive computer effort involved in obtaining the optimal inputs. 
Furthermore, the input matrix usually contains unknown parameters. 

In this case, Chen's optimum solution would not be of the bang-bang 
type and would be still more difficult to obtain. For all of these 
reasons it was concluded that at the present state of optimal input 
design methods an attempt to compare our selected inputs with an 
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"optimum input" would not be practical. Instead, a more limited 
approach has been taken described in the following section. 

4.4.2 Optimal Data Utilization for given Input Function 

We first point out the difference between the continuous and 
the discrete case. In the Maximum Likelihood (output error method 
for zero process noise) using the Newvon-Raphson approach with quasi- 
linearization, one obtains for the parameter update increment the 
following expressions: 


Continuous case: 



Discrete case: 


<1*) T 
'ae 'i 


(— ) 

^ae i 


fc’I 




) 


(115) 


The Cramer-Rao lower bound has been defined only for a vector of 
sampled measurements and not for the continuous case (2) and (18). 

For high sampling rate, one can define an approximate differential 
equation for M from equation. 115 in the following way: 


Set 6 ( 3v/3 0 ) ^ 


then 

N T -1 „ 1 

M*E S.R S i - 
i«l 


N i 

£ 5. R S At 

i = 1 1 


(116) 


( 117 ) 
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As N increases, At gets smaller and the right hand side of 
equation 117 can be approximated by 


M » (1/At) 



S T r" 1 S dt 


( 118 ) 


M 


-1 



(119) 


Taking the derivative of M - *- with respect to tf: 
dM^/dt* = -M" 1 OM/3t f )M‘ 1 

4 A 

or with equation 118 

dM* X /dt f = - ( 1/At ) M* 1 S T R* 1 S H -1 (120) 

The point is that even in a continuous formulation the time increment 
At between samplings must occur. Equation 119 is the correct formu- 
lation for the Cramer-Rao lower bound of the covariance matrix for the para- 
meters. (36) has a recursive formulation corresponding to equation 120. 

We can now use the approximately valid differential equation 120 
to obtain some insight into ways of best data utilization. Let us 
assume that we wish to prescribe certain values for the parameter 
standard deviations Oj and that we wish to compare the Cramer-Pao 
lower bound with these standard deviations. Since we are dealing 
not with the unknown actual parameter covariances but only with their 
lower bounds, we should apply some conservatism to the selected 
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that is we should select smaller than we really need for the 

specific data processing case. We thus require 

0 < Mtf (i, i) < o* (121) 

whereby is the value of M“1 at time t£. For non-zero 

values of S, the right hand side of equation 120 is negative 
definite and hence (i, i) are monotonically decreasing 

functions of t^. There will, thus, be a minimum time for which the 
constraints of equation 121 are satisfied. 

Another way of reducing the amount of measured data for the 
parameter identification is to select for the data processing those 
time periods for which the components of the matrix 

S T R~1 S 

have significant values. It follows from equation 119 that the 
Cramer- Rao lower bound then will ba particularly small. The 

components of also decrease with decreasing time element 

At between samples. 

Since it is impractical to use for the integration of equation 120 
infinity as initial condition, it is recommended to determine NT* 
for a small time period, say for N - 10, from equation 119 and 
integrate equation 120 with the solution to equation 119 as initial 
conditions. Since S includes parameter estimates, one needs a 
preliminary estimation of the unknown parameters in order to use 
equation 120. 
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4.4.3 Application to a Case of Lifting Rotor Parameter Identification 

The simplest mathematical model of the single blade flapping 
equation as expressed in equation 52 is used to identify the two 
unknown parameters: the collective pitch angle 0 O and the 

equivalent Lock number y. The angular acceleration, in the pitch 
stirring transient, is assumed to be J> = .1/ir and t = 12. 

Here we are concerned with the problem of designing the tests 
in such a way that the test data will be sufficient to determine 
the two unknown parameters y and 0 Q with good accuracy, i.e., 
to determine a suitable value of T that allows an accurate identi- 
fication of parameters. 

The simulated identification analysis was performed under the 
assumption of a random zero mean white noise sequence superimposed 
on the analytical flapping transient. This transient was obtained for 
0 O * 2°, y = 0.4 and y = 5.0. For convenience, the parameters 
6 = y 0 O and y instead of 0 O and y were identified. 


System and measurement equations corresponding to equations 103 
and 104 are: 


• 

r v 

x i 

1 

0 1 


V 


0 


S 






• 

L x 2j 

1 

-(<•>! 2 + (y/ 2) K(t)) -<y/2)C(t). 


- X 2- 


J Y/2)m 0o . 


( 122 ) 

♦ v(t) (123) 



where 

E (v(t)) e 0 E[v(t) v(t)] * 0.2 4(t - t) 


and [x x x ? ] * 


re 8] 
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We first sfcov in Table 11 the effect of data length on the 
parameters and their associated M~* (i, i) values. The iteration 
of the Maximum Likelihood method was begun with a 20 percent error 
in the parameter values. It is seen that a data length of t = 12 
- 14 is quite inadequate, a data length of t = 12 - 18 gives 
reasonably good parameters, while a data length of t = 12 - 24 
is much better and leads to a very small lower bounds of the parameter 
covariance matrix. Figure 14 shows the correct Happing response 
together with the simulated measurement data. Pitch stirring is 
initiated at t = 12. Figures 15 and 16 show M”*(y) and M“^(6) 
from equation 120 between t = 16 and t = 24. Two curves are plotted, 
one for the initial crude estimate of the parameters (y = 4, 6 = 8), 
and one for the final estimate of the parameters for t = 24, (y = 

4 . 91 , 6 = 9 . 83 ). The two curves are in this case not much different. 
Note the steep descent of the curves to about t = 17.5. It would, 
therefore, not be acceptable to use the data up to less than the 
time t = 17.5. However, there is another descent to t = 23.0, 
causing the improvement shown in Table 11. From Figures IS and 16 
it is clear that the selection of T = 24.0 is a good one, that the 
use of fewer data would result in substantial decrease in parameter 
accuracy, and that the use of additional data is unnecessary. 
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Figure IS. Plot of the Cramer-Rao Lower Bound of the Parameter covariance 
for the Parameter y from the continuous Formulation given by Equation 120. 
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Table 11. Parameter Identifiability for Different Data Length — 
Single Blade Model (Equations 122 and 123 ) with 
Parameters v and 6 A 0^ v . 

at 0 ' 



Initial Estimates 

4.00 

8.00 



Iteration 1 

4.29 

9.73 

48.0 

6.5 

2 

4.17 

9.71 

37.0 

8.1 

3 

4.10 

9.67 

37.0 

8.0 


ii) t = 12- 18 



Initial Estimates 

4.00 

8.00 



Iteration 1 

5.36 

9.67 

.096 

.032 

2 

5.23 

9.73 

.100 

.035 

3 

5.23 

9.73 

.094 

.035 

t = 12 - 24 
Parameter 

Y 

5 

M - 1 fy) 

M _1 (6) 

Initial Estimates 

4.00 

8.00 



Iteration 1 

4.94 

9.69 

.007 

.013 

2 

4.91 

9.85 

.008 

.015 


3 


4.91 


9.83 


.008 


.015 
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An analogous analysis was made for a forward flight condition, 
assuming the same pitch stirring transient and the same measurement 
noise. Now seven instead of three parameters must be identified. 

Figure 17 is a plot of the standard deviation (square roots 
of the Cramer- Rao lower bound) of the various parameter estimates 
versus the time of identification. The standard deviation is plotted 
as a percentage of the parameter value. The graphs are drawn under 
the following specifications: 

Parameter values chosen: 


0 j = .5 

M 


0 * 

M 

O 

O 

— 

v o 

0 6 

s 

®2 ® * 9 
0 j = - .8 

a. = x 

^ <£• 

Ml Ml 


0 e 2 e 3 
0 e 4 e 5 

V I 

V II 

** 

7 - s s, 

7.5C l 


e 5 * - 5 

0 6 * 1.0 

67 = Y = 4.9 

Pitch stirring excitation: 

0 j ■ 1.5 sin [a>(t - t 0 )} 
Ojj * 1.5 cos [«(t - t 0 )l 


(124) 


(125) 


where a 


0 

-( 0 . 1 )(t - t 0 i/» 


1 < *o 
t < t < T 


Measurement noise statistics: Mean * 0 

Std. deviation « .05 

Sampling time At » 0.1 time units 


Advance ratio v> « .4 
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The plot gives us the degree of identifiability of the parameter 
as a function of time length of identification. The Lock number y * 
e? is identifiable with a much higher degree of accuracy than the 
various parameters in the perturbation downwash equations. This was 
seen clearly in the simulation studies for the identification of the 
parameters. Beyond a time length of T * 18 the curves flatten out, 
indicating that measurement beyond that time does not improve the 
accuracy of the parameters identified. This is in agreement with the 
result for zero advance ratio. The above time length seems to be 
necessary and also adequate for identification purposes, for the 
given sampling rate and excitation. 

In our previous simulation studies we have used a sampling 
time of T = 12 time units. From Figure 17 it appears that 
inaccuracies in our identified parameters could be attributed to 
inadequate data length for identification purposes. This factor 
will be taken into consideration in parameter identifications using 


test data. 
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5. BRIEF DESCRIPTION OF THF. EXPERIMENTAL SET-UP 


The basic purpose of the experimental set-up is to measure 
continuously the flapping response of the rotor blades to cyclic 
pitch stirring excitation. A sampled length of the response, together 
with the input excitation, is used to identify the dynamic inflow 
parameters and hence determine the feedback effect of the rotor 
wake on the flapping response. 

The experiment can be split up into three independent circuits: 

1. The strain gauge (the flapping response measurement) circuit 

2. The pitch resolver circuit 

3. The rotor resolver circuit 

A brief description of the above circuits are given below. 

For a detailed description of the experimental equipment and 
procedure see (12). 

5.1 STRAIN GAUGE CIRCUIT 

Four strain gauges are mounted on the flexure of each blade 
to form a Wheatstone bridge. Two strain gauges on the top and 
two strain gauges at the bottom of each flexure are so connected 
that the torsional and the lead-lag motions, if any, are annulled. 

The rotor is considered to be very stiff in torsion and lead- lag. 

A schematic diagram of the strain gauge circuit is shown in Figure 18. 

Power is supplied to two arms of the bridge through two slip 
rings. The signal is taken out through two other slip rings from the 
other two arms of the bridge. The signal is passed through medium- 
gain amplifiers and recorded on a six channel FM tape-recorder. 
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Figure 18. Schematic Diagram of the Balancing Network for each 
Strain Gauge Circuit (S.G. - Strain Gauge; S.R. - Slip Ring). 
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5.2 PITCH RESOLVER CIRCUIT 

The resolver is a rotating transformer whose output varies sinu- 
soidally • with the angular position of the resolver shaft. The 
resolver shaft is connected by means of a sprocket drive to the inner 
shaft Figure 19 . The relative velocity of the rotor and the inner 
shafts provide the pitch stirring excitation. A brief description 
of the drive mechanisms will be given later. The input to the 
resolver is an oscillator whose frequency and amplitude can be varied 
to get desirable output signal from the resolver. The output signal 
is of a varying amplitude depending on the angular position of the 
resolver shaft with a carrier frequency corresponding to the 
input (oscillator) frequency. The signal is then passed through 
a full-wave rectifier-low pass filter circuit to remove the carrier 
frequency. The signal is then fed into a low-gain amplifier to adjust 
the output level of the signal. It then goes through a 3-way, 2-position 
switch which works as follows (Figure 20): 

When the switch is in the "OFF" position, the input to the 
recorder is a 2.5 volt D.C. battery signal. When the switch is 
flipped to the "ON" position, the following events occur simultaneously: 

a. The motor drive to the inner shaft is activated 

b. The signal sent to the resolver is no.- the resolver signal 

c. The solenoid that retains the inner shaft at a fixed (trim) 
condition, is released 

5.3 ROTOR RESOLVER CIRCUIT 

An oscillator provides the input to the resolver which provides 
the angular position of the main rotor shaft in a manner similar to the 
input to the resolver on the inner shaft . The output of the rotor 
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Figure 20. The Resolver Circuit of the Pitch Stivring Excitation. 
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resolver is then passed through a full wave rectifier circuit. 

It is then input into a low gain amplifier for output level adjustment 
as before. It is then sent to a two-way, two-position switch, the 
output of which is input to the tape recorder. In the "OFF" 
position the tape recorder input is a 2.5 volt D.C. signal. When 
the switch is flipped to the "ON” position, the motor that drives 
the rotor shaft is energized and also the rotor resolver signal 
is input to the tape recorder. A schematic diagram of the circuit is 
shown in Figure 21. 

5.4 PITCH STIRRING EXCITATION 

The inner shaft is a cylindrical rod which passes through the 
hollow shaft with an eccentric pin mounted at the end as shown in 
Figure 19. TWo sets of pitch control flexures are mounted on the pin. 
Each set of flexures is clamped to two opposite blades. Rotation 
of the inner shaft is effectively the same as rotating a tilted 
swash plate. The collective pitch of the blades is adjusted by 
loosening the pitch lock nut adjusting the blade pitch on the pitch 
screw and then relocking the pitch lock nut. The drive is mounted 
on the other end of the inner shaft. Basically two drive mechanisms 
are used: 

1. A motor is used to drive the inner shaft using a sprocket 
drive. The acceleration of the motor provides the required 
transient excitation. This is seer in Figure 22a. 

2. A coil spring mounted on the base of the inner shaft is 
used to drive it through 90, 180 and 270 degrees. The 
acceleration of the shaft would be proportional to the 
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Figure 21. The Rotor Resolver Circuit. 
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amount that the coil spring was wound prior to its release. 

A view of the spring excitation is shown in Figure 22 b. 

5.5 BRIEF DESCRIPTION OF THE EXPERIMENTAL PROCEDURE 

Ihe stands are adjusted so that the rotor plane is parallel 
to the base of the wind tunnel for study of the flapping response 
for zero angle of attack. The collective pitch of each blade is set 
in the following maimer: 

In order to make relative collective pitch changes, a beam 
of light is focused on the small mirror glued to the root of each 
rotor blade as shown in Figures 23a and 23b. At zero collective 
pitch setting, the beam striking the mirror at an angle of incidence 
of 30° is reflected to a position marked X* on the calibrated 
scale. On changing the collective pitch Figure 23b, the angle of 
incidence of the beam changes, which is reflected to the position 
marked Xj, which directly reads the change in the collective pitch 
setting from the position X*. 

To set the collective pitch to zero degrees, the zero degree 
eccentric is mounted on the inner shaft. The pitch angles of the blades 
are adjusted till each of the blades have a minimum flapping response. 
This is studied on the scope. 

A qualitative judgement regarding the relative accuracy of the 
collective pitch settings is made by using the stroboscope. A 
photocell reflecting off thin reflecting strips (corresponding to 
each blade) on the rotor shaft is used to trigger the stroboscope. 

The blades can be observed to have the same flapping angle. If not, 
small adjustments have to be made in the collective pitch settings. 
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The zero eccentric on the inner shaft is replaced by the 
± 1.4 degree eccentric. The inner and the rotor shaft resolvers 
are now synchronized. This is done by first locking the rotor shaft 
and the inner shaft by a locking pin. A screw on the rotor shaft 
which is located at the zero azimuth angle is used to activate 
a magnetic pick-up which generates a voltage blip. The inner 
and the rotor shafts are run in the locked position. The rotor 
and the inner shaft resolvers are adjusted such that their zeros 
pass through the center of the blip of the magnetic pick-up. The 
inner and the rotor shaft resolver positions are now synchronized. 
Since the resolver output amplitudes are not exactly sinusoidal for 
constant rotational speed, the resolver signals of both the inner 
and the rotor shafts are recorded with the two shafts locked. This 
provides information for calibration of the resolvers. 

The shafts are uncoupled and the drive for the inner shaft is 
set up. The strain gauge circuits are balanced as described before. 
The flapping deflection has now to be calibrated. This is done by 
using the following blade specifications: 

(a) one inch of tip deflection corresponds to 6.23° flap 

(b) one inch of tip deflection requires a moment of .562 lbf-in 

With the above information, a 10 gram mass at 7.17 inch (blade 

tip) distance is found to create a deflection of 1.746° flap. The 
amplifier gain is adjusted for an output of 1.746 volts, thus giving 
a one volt/degree flap deflection. The tape recorder is zeroed 
and then calibrated by recording two known levels of D.C. voltages. 
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The rotor structure is mounted on a swiveling base so that the 
rotor shaft can be tilted for different angle of attack positions. 

The rotor shaft tilted at « = -3° is shown in Figure 24a and 24b. 

A fairing is mounted on the rotor shaft for a streamline flow past it 
and to avoid the effect of the rotor shaft on the blade wake effects. 

The motor that drives the rotor shaft is turned on by the switch 
MS-2 (Figure 21). The inner shaft is held in position mechanically 
by means of a solenoid. The rotor shaft speed is adjusted by a 
potentiometer, thereby changing the rotational flap-bending stiffness 
u»j of the rotor blade. The wind tunnel speed is adjusted for tne 
required advance ratio. The tunnel speed is measured by a manometer. 

A six channel FM tape recorder, which measures the response of the 
four blades and the two resolver signals, is turned on to record the 
signals. Switch PS-3 (Figure 20) is in the closed position and PS-2, 
the inner shaft drive motor, is in the "OFF" position. When the 
switch PS-2 is flipped, the inner shaft drive motor is turned on and 
the solenoid that holds the inner shaft is released at the same time. 
The acceleration of the inner shaft drive provides the transient 
excitation. Switches PS-3 and PS-2 are turned off in sequence and 
the recording stopped. A view of the experimental set-up is seen 
in Figui'e 25. 
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6. EXPERIMENTAL DATA ANALYSIS 

Hie transient test results together with the trim conditions 
are recorded in analog form on magnetic tape. Hie two resolver 
signals and the signals from the blade flapping strain gauges are 
processed on a PDP-8 and PDP-12 mini-computer complex. Prior to the 
identification analysis, several preliminary x*aw data manipulations 
are necessary. Hie analog data are digitized, the data is transformed 
from rotating to fixed coordinates and the transient perturbation 
results are separated from the trim values. Hie small bias values in the 
blade flapping measurements are removed by identifying the bias values 
in each identification run. 

For low advance ratios (y{.4), it is seen (from 4.3) that the 
multiblade formulation gives accurate response with the periodic 
terms left out. When the periodic terms are neglected, the 0^ 
equation 67 becomes uncoupled from the rest of the equations. Hence 
the response can be neglected from the identification procedure. 

Even for m = 0.6, a comparison of the y* identification between 
a model with the periodic terms, the 0^ equation and reverse flow 
effects included and a model neglecting all of the above effects 
is surprisingly accurate . These are shown in Table 12 

and Table 13. Thus it is seen that neglecting the reverse flow and 
the periodic terms at low advance ratios is quite acceptable. At 
low advance ratios the reverse flow regions are restricted near the 
rotor hub, thus having negligible effect on the flapping response 


of the rotor blades. 
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Table 12. - 13. Comparison of the Identified y* Values Using a 
Mathematical Model which Includes Reverse Flow, 
Periodic Terms and 8 d equation (Table 13) with 
a Model Excluding the Above Effects (Table 12) 
Using Data at y * 0.6. 


Table 12 



Parameter 

* 

Y 

Least 

Squares fit (R£ 
o 6j o &n 

Initial 






Value 


3.00 

.01735 

.39835 

.06648 

Iteration 

1 

3.529 

.00184 

.00812 

.00767 

2 

3.654 

.00163 

.00859 

.00754 

3 

3.654 

.00163 

.00859 

.00754 

4 

3.654 



Table 13 



Parameter 

* 

Y 

Least Squares fit (R) 

2 2 2 
o 8 0 a Bj o 8 jj 

2 

° e d 

Initial 







Value 


3.654 

.00112 

.00816 

.00716 

.01540 

Iteration 

1 

3.630 

.00112 

.00866 

.00719 

.01542 

2 

3.595 

.00112 

.00795 

.00726 

.01557 

3 

3.584 

.00112 

.00781 

.00725 

.01562 


3.580 
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An interesting observation is that at higher advance ratios the 
bias terms are rather large. This is because at high y, the 
trim flapping response has (seen experimentally) a small variation 
of the amplitude of the trim condition with a large time period. This 
can very well cause the trim subtracted from the transient response 
to be slightly different from the actual trim at which the transient 
was measured, thus giving rise to the bias values. 

If all the blades were identically set, the trim values for 
the flapping response of the blades in the non-rotating system would 
be nearly constant . From equation 67 it is apparent that 8^ has 
a moderate 2/rev. in its trim condition, 3j and have a smaller 

4/rev. input effect, and 3 q has a very small 4/rev. for its trim 
response. A typical average trim condition data is given in 
Figure 26. 

If a constant approximation is used for the trim values, it 
will correspond to the starting values of the transient responses. 

This would ensure response bias values of approximately zero, though 
the 2/rev. and the 4/rev. trim conditions show themselves in the 
transient data. If the periodic trim conditions for data obtained 
at low advance ratios are used, then the transient responses will be 
rid of both the bias and the 2/rev. and the 4/rev. trim variations. 

6.1 PITCH STIRRING EXCITATION 

Two rates of acceleration of the inner shaft are obtained by 
adjusting the potentiometer setting of the eddy current to the motor 
that drives the inner shaft. Plots of the slow and the fast excitation 
are shown for normalized 9 jj as a function of non-dimensionalized 
time in Figure 27 and Figure 28 respectively. Care was taken to use 
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as a Function of Non-Dimenslonalized Tims. 



125 




126 


a length of data in which the excitation velocity is uniformly 
increasing. For both rates of excitation, the frequency of excitation 
becomes approximately constant beyond t * 30 and any data utilized 
in this range would tend to give biased parameters whose values 
depend on this constant excitation frequency. 

The inputs are mostly progressive pitch stirring accelerations 
of the inner shaft. Data is also gathered using regressive pitch 
stirring acceleration and spring loaded pitch stirring excitation. 
Details of these inputs are given in (12) . These data sets are used 
for verification of the identified model response prediction. 

6.2 y* IDENTIFICATION RESULTS 

Figure 29 is a plot of y* versus advance ratio for different 
values of the collective pitch setting at 1.18. The values of 

y* from momentum theory (6^0°, equation 58) are plotted for comparison. 
This analytical result is for low collective pitch settings. As the 
advance ratio increases, the rotor wake gets washed away faster and 
hence, the feedback effect of the dynamic inflow on the flapping 
response decreases. This is seen in the increasing value of the 
identified y* with advance ratio for all collective pitch settings. 

Figure 30 shows similar results for = 1.24. The analytical 
model approaches the true value of y asymptotically. The experi- 
mental results for 0 O * 0° shows that around p » 0.75, y* is 
approximately equal to the value of the blade Lock number y. This 
shows that at this high advance ratio, the effect of the rotor down- 
wash has become negligible. 
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Figure 29 and Figure 30 are combined and plotted in Figures 31a 

and 31b, where y* is drawn as a function of Q q , the collective 

pitch setting. The consistent trend seen in these curves is that, 

with increasing value of 0 O , the y* value first drops and then 

increases beyond a collective pitch setting of around 3°. This is 

in apparent contradiction to the equivalent Lock number formulation 

given by equation 57, which indicates that, since with increasing 

collective pitch 0 Q the induced mean downwash v monotonically 

increases (11), the equivalent Lock number increases. The discrepancy, 

however, is due to the fact that the tip loss factor B Is decreasing 

with increasing 9 q , causing an apparent decrease In the value of the 

d 

Lock number which varies as B . 

* 

Also plotted in Figure 29 and Figure 30 are graphs of y for 

trim conditions of 0 = 3.5°: « = -3°. This increases the downwash 

o 

over the trim condition of 0 Q = 3.5°; <* = 0°, and hence, decreases 

the effective angle of attack of the blade. This is seen in the 
higher value of y*. This agrees with the trend of y*versus 0 O 
as found before. 

Studies have also been done to determine the effect of a 
range of a from 0° to -6°. The shaft tilt corresponds to a combination 
of collective and lateral cyclic pitch change of the rotor blades. 

This is because the changes in the angle of attack (AO) due to a 
shaft angle of attack of o ar«: 

At p = 0 and ir (fore and aft positions) : 

A0 = mo/ ( r/R) (126) 

At p = n/2 (advancing side) : 

A0 = Mo/(r/R+y) (127) 

At p«3n/2 (retreating side) : 


A0 = ya/(r/R-M) 


(128) 
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Figures 32a and 32b are plots of Ay* versus shaft tilt angles. 
Since these are essentially plots of y* versus collective pitch 
change (on a different scale) , the trend of these graphs should 
agree qualitatively with those in Figures 31a and 31b. This is an 
independent verification of the results shown in Figures 31a and 31b. 

y* has been found to have very good identifiability. Data 
length study has been done for transient data at wj = 1.24, 6 0 = 5°, 
y * 0.3. y* was identified using seven transient data lengths ranging 
from two to five rotor revolutions. In all of the cases, the y* 
identified to within 4% of one another. 

6.3 L-MATRIX MODEL IDENTIFICATION 

A simplification in the flapping response of the L-matrix is 
made when the g Q equation is neglected from the set of equations 
used for identification (see equation 64) . This is done for the 
following reason: 

Numerous identification tests have shown that both and 

t q A 1/Kjj have poor identifiability (see equation 71). These are 
the two parameters associated with A Q . The primary effect of 
A 0 is on the g Q response and the g 0 transient response is small in 
comparison to gj and gjj responses. Since g Q is weakly coupled 
with gj and g j j responses at low advance ratios, the equations 
governing the g ft and the \ Q responses can be neglected. 
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The reason for the small 0 O transient response, and consequently 
poor identifi ability of the two associated parameters Ljj and t 0 seems 
to be due to the lack of collective pitch excitation. Pitch stirring 
excitations in the pitch and roll directions generate adequate 0 j and 
*11 responses to identify the parameters associated with them. 

The above simplification reduces the identification problem to 
one of identifying the four elements L 22 » L 33 . L 32 and L 23 of the 
L-matrix, xj, xjj (see equation 71) and y, based on the measure- 
ments and Bjj, 0j and 0jj. y appears as a product with B 2 , 
and B* in the flapping equations. The tip loss factor B is not 
known accurately as a function of the blade pitch angle. Simulation 
studies have shown that y identifies very accurately (within 3% 
of its true value). Hence the value of B is assumed to be 1.0 and 
the Lock number y is assumed as a parameter to be identified. The 
identified y will thus represent, not the true value of the Lock 
number, but a product of the Lock number and a power of 0 between 
2 and 4. Roll and pitch time constants Tj and xjj are assumed 
to be identical in theory ( 11 ) , but with increasing advance ratio they 
are expected to have increasingly different values. 

Figure 33 is a plot of y versus 0 Q . The trend of the curve 
With increasing 0 o , according to steady momentum theory, should be 
due to decreasing value ef the tip loss factor. 
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But the value of the Lock number increases and then drops beyond a 
certain 0 O , This trend is also found in a similar plot for hovering 
analysis (12). The plot also shows a decreasing variation of y with 
increasing advance ratio (graph drawn for y ■ 0.2 as compared to the 
one drawn at y » 0 . 1 ). 

Figure 34 shows a plot of the parameter L 22 as a function of 
the advance ratio for different identification runs. The corresponding 
values predicted by momentum theory are plotted for comparison. 

The experimental values are several times larger than those from 
momentum theory. is similarly plotted for two different otj 

values in Figure 35. 

The parameters L 2 3 and show no consistent trend with 

variation of the trim conditions. The values from all the experimental 
results lie around the value of zero, as predicted by momentum theory. 
All identified values of L 2 3 are between -0.15 and +0.15. 

All identified values of L 32 lie between -0.15 and +0.4. 

From several studies, it is seen that Tj and t jt have approxi- 
mately equal values at y = 0.1 regardless of their trim conditions. 

The value is very close to the theoretical value in (9 ). At y « 0.2, 
the value Tj becomes larger and tjj becomes correspondingly smaller. 
The ratio of tj to Tjj ranges from 1.5 to approximately 2.5. A 
typical conparison is shown in Table 14. 
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Figure 34. Plot of Identified Parameter L 22 Versus Advance Ratio. 
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6.4 COMPARISON OF IDENTIFICATION AND PREDICTION STUDIES 

As explained before (page 126), the transient data length used 
for identification is selected such that the input transient has 
a uniform increase in velocity. Progressive input excitations of 
the slower type (Figure 27) and their corresponding responses are 
generally used for the identification of the parameters. 

Once the parameters have been identified, based on a certain 
length of data, the goodness of these parameter values, and hence, 
that of the corresponding mathematical model, has to be ascertained. 
Prediction studies are, hence, required, and are made in the following 
manner: 

The identified parameter values are -nserted in..; the mathematical 
model. The model is now complete. This complete model together 
with the faster transient input (Figure 28) is used to determine 
the response of the mathematical model to be compared with the response 
of the experimental model to the same faster excitation. This 
prediction study is done for closely similar trim conditions as 
those in the corresponding identification study. 

Typical examples of studies done are given here. 

1. In Figures 36a, 36b and 36c, data set with 9 q = 0°, oij = 1.24 
and y = 0.4 is modelled using the y* model and the L-matrix 
model programs. This data is modelled so well with the y* model 
that there is hardly any improvement by using the more complete 
L-matrix model. 

2. Prediction curves for data with u, = 1.116, 9 =5° and u = 0.1. 

1 o 

A comparison of the three models is studied for prediction 
results in Figures 37a and 37b. The three models are: 
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(a) y model (neglecting downwash) 

(b) y* model (obtained from data with o>j ■ 1.24, 0 q = 5° 
and n ■ .1) 

(c) the more complete L-matrix model (parameters obtained from 
data with o>j * 1.24, 0 Q * 5° and p * 0.1) 

Though the model without downwash shows high r. m. s. fit errors, 
the ether two models give good degree of fit prediction. 

The y model gives larger amplitude responses as compared to the 
other two. Hie Lock number y is the ratio of the aerodynamic 
forces to the gyroscopic (inertia) forces. Increase in the aero- 
dynamic forces increaseethe amplitude of the response.' Increase 
in the inertia forces decreases the amplitude of the response. 

Hence, with increasing value of y, the response has an increasing 
amplitude. This is also verified by the flapping equation 52 
which shows that, since y multiplies the forcing function, the 
response is directly dependent on it, even though the damping of 
the system is increased by increasing the value of y. 

Frequency response curves show that the y* model has a decreasing 
value of the anqalitude response with increasing progressive 
frequency excitation (beyond a certain frequency) , whereas the 
true response has an opposite trend (37) . This gives rise to 
greater discrepancies at higher progressive frequencies between 
Y model and the physical model. This is clearly seen in the 
$2 and 0jj prediction curves, since the prediction studies are 

made at a much higher frequency of input transient excitation as 

* 

compared to the frequency range in which the y model was determined. 
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The initial conditions for these prediction studies were chosen 
at their corresponding identified values. 

3. Prediction curves for data with wj * 1.18, 9 0 * 5° and p « 0.2. 

The L-matrix parameters and y were chosen from the data set with 

(i>j ° 1.24, 9 q ■ 5° and p » 0.2. The three models are studied as 

before (see Figures 38a and 38b)* All the features found in the 

previous prediction study with p * 0.1 are found in this study. 

An important observation is that the model neglecting the dynamic 
* 

inflow (y model) has a significantly smaller error at p 0.2 
as compared to that at p * 0.1. This indicates the fact that 
with increasing advance ratio, the feedback effect of the down- 
wash is diminishing. 

4. In Figures 39a, 39b and 39c, data set with 9 Q = 5°, Wj = 1.24 

* 

and p - 0.2 is modelled using the y model and the more complete 
L-matrix model (with diagonal terms) . The plot shows a signifi- 
cant improvement in the fit by using the complete L-matrix model. 
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7. CONCLUSION AND DISCUSSION OF RESULTS 

Based on the preceding analysis of the feedback effects of the 
rotor dynamic inflow on the flapping response of the rotor blades to 
accelerated pitch stirring excitation, several interesting results follow. 

7.1 CONCLUSION OF METHOD ANALYSIS 

Apart from the several theoretically justifiable properties of 
the Maximum Likelihood method as outlined in section 2.7, simulation 
studies and comparison with other identification techniques give the 
following results: 

1. The Maximum Likelihood method works well for both the single 
blade and the multiblade application in simulation studies. 

2. The Cramer-Rao lower bounds for the parameter covariances 
obtained from the Maximum Likelihood analysis provide a good 
measure of accuracy of the identified parameters and are clearly 
superior and more meaningful than the covariance estimates 
determined with other known methods. 

7.2 RESULTS OF THE SIMULATION STUDIES 

* 

1. Single blade identification of y and 0 q together with the Initial 
values of the flapping response at p B 0.4 give good results. This 
analysis is not used for experimental data because of small 
differences between the different blade responses. 

2. The parameters of the entire L-matrix (as given by equation 96) 
converge in the simulation identification of the perturbation 
flapping response model, though with limited accuracy. The 
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off-diagonal terms and L^ end the time constants t q , tj 

and tjj have been found to provide poor identifiability in certain 

cases. 

a 

3. . The Y Identification study provides consistent and fast convergence 
under all test conditions. 

4. Eigenvalue analysis of different forward flight models show 
(Table 7) that for the cases studied, at low advance ratios (y<. 4), 
the error in neglecting the periodic terms in the flapping equations 
is negligible (< 2% in the eigenvalue comparisons). On the other 
hand, neglecting the inflow in the mathematical model gives rise 

to significant errors. 

5. The eigenvalue variation with advance ratio (Table 9) is negligible. 
This is based op no variation of the inflow parameter values 

with advance ratio. 

6l. At high, advance ratio (u»0.8), there is negligible change in 

the eigenvalues (Table 10) vhe$ the ffftdb^ck effects of the inflow 
is neglected. 

7.3 RESULTS FROM EXPERIMENTAL DATA ANALYSIS 

1. The cyclic pitch stirring excitation with approximately constant 
stirring acceleration is adequate for identification of the 

it 

parameter y (equation 72) and some of the L-matrix parameters 
(equation 74). To identify the parameters associated with the 
equation accurately (equation 74) , collective pitch excitation 
is probabl^ required. 

2* Progressive transient excitation data gives identifiable parameters* 
The same identification study on regressive cyclic pitch stirring 
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transient excitation data gives equally good results. 

* 

3. y is an accurately identifiable parameter for all trim conditions 
(the collective pitch setting, the blade rotational stiffness, 
the rotor shaft angle of attack and advance ratio). This concurs 
with the findings of the simulation studies. In contrast, the 
L-matrix parameters identify only upto an advance ratio of 
M s 0.4. Beyond this advance ratio it is not possible to identify 
these parameters. 

4. The identifiability of the parameters from the experimental data 

is in dose agreement with the findings in the simulation studies. 

The time constants t q , and t ^ have increasing identifiability 

in that order. Their effects on the flapping responses are 

lower than those of the parameter!, of * ' L-matrix. 

* 

5. The y model adequately represents the feedback effects of the 

induced inflow on the flapping response for advance ratio p*0.4. 

At low advance ratios (u<0.4), the L-matrix model provides a 

much better correlation to the experimental flapping responses 

* 

as compared to the y model. With increasing advance ratio, the 
difference between the fits of the two mathematical models and 
the experimental results decreases. 

6. The identified values of some of the L-matrix parameters agree 
reasonably well with those obtained by using momentum theory 
(Figures 34 and 35). Most of the Identified values of L 22 and 
Ljj fall between the values obtained from momentum theory and 
the empirical values given in (10). Empirical L 22 and values 
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are approximately -0.5 and -0.2 respectively for :dvence 
ratios between 0 and 0.4. 

7. and Identify to values around the momentum theory value 
of zero. They show no consistency or trend in their variation 
with the trim conditions. The identified values of lie 
between -0.15 and 0.15; lies between -0.15 and 0.4. 

8. The accuracy of the identified mathematical model is determined 
by applying it to tests not used for the identification. Examples 
shown in Figures 37 a-b and 32 a-b show the adequacy and goodness 
of the inflow models chac have been used. 

7.4 AREAS OF FURTHER STUDY 

1. The plots of Figures 31 a-b show a consistent trend that has not 

been explained T»y theoretical models (see 6.2). This suggests 

the need for a more complete mathematical model to represent the 

* 

variation of the equivalent Lock number y with trim condition. 

2. The parameters associated with the perturbation equations for 

v q (equation 74) are not identifiable using cyclic pitch stirring 
excitation. In order to make a complete identification, collective 
pitch transients will have to be added. 

3. In the preceding study, dynamic rotor inflow is coupled to blade 
flapping only. It is unlikely that blade lead-lag or torsional 
flexibility will have a substantial effect on the dynamic inflow. 
However, the effect of the dynamic inflow on lead-lag and torsional 
deflections is expected to be of importance and should be studied. 
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8.1 NOMENCLATURE 


B 

B(j/k) 

C(t) 

C T 

Si 

S. 

D(t) 

E 

E(0/Z) 

F(t) 

G(t) 

H(t) 

I 

J 

Si* v % 

K(t) 

K(j> 

L 

L ll* L 12’**’ 
M 

Sc 

P(t) 


blade tip loss factor 

innovation covariance matrix at time t, given 
measurements till time t^ (equation 32j 

aerodynamic damping (equation 33) 

rotor thrust coefficient, positive up (equation 130) 

rotor pitching moment coefficient, positive nose 
up (equation 130) 

rotor rolling moment coefficient, positive to 
right (equation 130) 

measurement input matrix 

energy term (equation 105) 

expected value for the probability density of 0 
given the measurements Z 

rotor state matrix 

rotor input matrix 

measurement state matrix 

identity matrix 

scalar cost criterion 

nondimensiona* apparent mass and inertia of 
impermiable disk 

aerodynamic stiffness (equation 54) 

Kalman filter gain at time 
rotor induced flow gain matrix 
parameters of the L * matrix (equation 73) 
information matrix 

moment at the rotor hub of the k th blade 

the state covariance matrix 

when K =K T =K T 
1 l I 2 1 
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NOMENCLATURE (continued) 




Q 

R 


T k 

U 

P 

U T 

W 

Y J 

z 

a 

b, N 
c 

ft.) 

h(.) 

g 

m Q (t), tt x (t) 

PC.) 


V T 


8 

t 


combined covariance matrix, where the parameters 0 
are included In an augmented state 

system equation noise covariance matrix 

measurement equation noise covariance matrix and 
rotor radius 

sensitivity matrix at time (equation 116) 
thrust at the rotor hub of the k th blade 
relative blade normal velocity (equation 133) 
relative blade tangential velocity (equation 132) 
positive definite weighting matrix 

set of observation vectors till time t^ (y^, y^ , • • » y^) 
set of all observation vectors (y^, y^, , y^) 
blade section lift slope 
number of blades on the rotor 
blade chord 

function of variables in parenthesis in the 
system equation 

function of variables in paranthesis in the 
measurement equation 

vector of measurement bias 

aerodynamic flapping moments (equations 55 and 56) 
probability density function 
initial time 
final time 

complex variable in the laplace transform 

nondimen8ional time for which the period of one 
rotor revolution is 2ir 


u(t) 


input control vector 
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NOMENCLATURE (continued) 


v(t) 

w(t) 

x(t) 

y(t) 

a 

®o» e d» e i» e il 


6 «Y© 

1 c 
«(t) 

A 


\>» a i# A 


II 


IC 

n(x) 

(a 

w l 

a 

v(j) 


V II 


v k 

P 

Y 

ra) 


measurement noise vector of covariance R 

system noise vector of covariance Q 

rotor state vector 

rotor output or measurement vector 

hub pitch angle, positive nose up 

flapping angle of the k th blade, positive up 

multiblade flapping coordinates: coning, differential 
coning (only for even-bladed rotors), longitudinal 
and lateral cyclic\ flapping . _ % . ' i 

arimuth angle of the k th blade 

parameter identified in the single blade analysis 

delta function 

nondimensional normal inflow 
uniform, longitudinal and lateral inflow components 
adaptation factor for the first mode representation 
first rotating mode shape of a rotor blade 
angular pitch stl.ring speed 

blade flapping natural frequency in the rotating system 
rotor angular speed 
rotor advance ratio 

innovation vector at time t^ (equation 23) 

uniform, longitudinal and lateral perturbation 
induced inflow components 

induced inflow of the k th blade 

air density 

blade Lock number 

system noise matrix (equation 3) 
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NOMENCLATURE (continued) 


V V e n 


® 1 * ®6 


0 

OF 


A 

C 

Superscripts 


T 

* 


collective* nose down cyclic end left cyclic 
pitch angles respectively 

instantaneous blade pitch angle 

parameters of the inflow perturbation model 
(equations 102 and 104) 

vector of unknown parameters 

rotor solidity ratio (bc/vR) 

standard deviation of the i th parameter 

positive definite weighting matrix 

measured states in the system equation 


time derivative 
estimated value 

mode of the probability density function 
transposed matrix 
equivalent value 
mean or trim value 


Subscripts 

m 


i 

k 


o» d. I, II 
E 


Symbols 


A 


measured variable 
i th sample of the variable 
blade number or iteration number 
multiblade variables 
empirical value 

approximate equality 
defined by 



162 


NOMENCLATURE (continued) 


I I 

<\» 

n 

A 

Z 

.a 


{ 

< 

| ug - u T wu 


determinant of a matrix or absolute value of a number 
equivalent to 

product of the given terms 
Increment 

summation of the given terms 

definite Integral with limits from a to b 

less than or equal to 

greater them or equal to 

euoledlan norm 


W 
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8.2 DERIVATION OF THE THRUST AND MOMENT COEFFICIENTS C T , C. 

A L 

Cjj AS A FUNCTION OF THE STATE VARIABLES 
The state of the ayetem Is 


* T ‘ ». K «I «I »II «II v„ Vj v„ 6 d 6 d l 



Cj^ ■ - (l/pirt^R 5 )^^ M^ cos 
k-1 

N 

Cjj • - (l/pirfl 2 R 5 ) sin * k 

k-1 

N 

Gj, - (l/pirft 2 R*) T k 

k-1 


AND 


(129) 


(130) 


We now determine the thrust end moment of each blade to be substituted 
Into the equations 130. 
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Blade pitch angla:* 9 ■ 0 - 6 T sin \p + 6^ cos 

(131) 

Induced flow: v ■ v Q + x cos ip + Vjj x sin i|> 

Relative tangential velocity: 

U T - x + w sin \p (132) 

Relative normal velocity: (4- are up) 

Up - v - x$ - uB cos * (133) 

The thrust and moment exerted by each blade, as given under the 
assumption of no reversed flow are: 


T k“ 


pacfl 2 R 3 

2 


/' 

•'o 


(U^ 9 + Up U T ) k dx 


(134) 


*k- 


pacfl 2 R t> 


/' 

J o 


(U^ 9 + Up U T ^x dx 


(135) 


Equation 132 can be substituted into equation 134 with the aid of the 
transformation 


3 k - 6 0 + B d (-D fc + 6 X cos i|> k + Bjj sin <|> k 

B k « & Q + B d ("l) k + + Bjjlcos <» k + (B n - Bjlsin «|» k 


( 136 ) 


The result, after some algebraic manipulation, is: 
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T fe - P - a £ jj 2R - {(B 3 /3 + By 2 /2 + B 2 Msin ^ - By 2 cos 2^ k /2)e Q 

- (B 2 y/2 + (B 3 /3 + 3B y 2 /4)sin * k - B 2 U cos 2^/2 

- By 2 sin 3^/4^ + ((B 3 /3 + By 2 /4)cos <J> k + B 2 y sin 2* fc /2 

- By 2 C 08 3^/4) 0 n + (B 2 /2 + By sin i|» k )v o 

+ (B 3 cos tJ> fe /3 + B 2 y sin 2 ^ k /4)v I + (B 3 sin * k /3 + B 2 y/4 

- B 2 ycos 2* k / 4)v II " ( ® 3/3 + B2w8in V 2)(B o + h (_1)k) 

- (B 3 cos * k /3 + B 2 y sin 2 iJ» k /4) (Bj + B n ) 

- (B 3 sin tl> k /3 + B 2 y/4 - B 2 ycos 2^/4) (B n - Bj) 

- (B 2 y cos * k /2 + By 2 sin 2^/2) (B Q + B d (-l) k ) 

- (B 2 y/4 + By sin 4» k /4 + B 2 y cos 2ty k /4 + By sin 3i|» k /4)Bj 

2 

- (By cos t[» k /4 + B y sin 2^/4 - By cos 3<|> k /4)0jj 

- < ((B 0 + B d (-l) k + Bj cos ip k + B n sin i|i k )(.4411 y 2 

sin 2i/> k + 1.0676 y cos t|> k ) + (B Q + B d (-l) k (137) 

+ (B x + B II )cos * k + (Bjj - Bj)sin ^ k > (y sin <|» k (-.2119) - .0838))} 
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Equation 130 £or C^, is now substituted in the previous equation 137. 
The signs of the Induced downwash v q , and v are changed for 
the proper positive down) convention. 

C j - (ao/2){(B 3 /3 + Bu 2 /2)e Q - (bVz^ 

- (B 2 /2)v q - (B 2 p/4)v n - (B 3 /i)B o b (By 2 /2)sin 2t 0 d 

- <B 2 m/4)0 u + k ((.4411)p 2 sin 2t 0 d 

-.5338y Sj. + -0838 B q + (.1059)y (E - Bj.))} (138) 

Equation 133can now be substituted in equation 135, and using the 
transformation equation 136and equation 130 » the following expansions 
for and C^, in terms of the state variables, are obtained: 

= -(l/pffi2 2 R 5 )(pacn 2 R 4 /2){((B 4 /4 + B 2 y 2 /4)0 Q 

- (B 3 U /3)0 i + (B 3 /3)v q + (B 3 y/6)v n - (3^/6^ 

b 

- (B 4 /4)(B q + B d (-l) k )) Y sin \ + «2B 3 y/3)0 o 

k=l 

- (B 4 /4 + 3B 2 y 2 /8)6 I + (B 2 y/2)v Q + (B 4 /4)v IT - (B 3 y/3)(B 0 + B d (-l)S 

b 

- (B 4, 4)(B - 6 r ) - (B 2 y ? 78)B T )(l/2> Y U-cos 2*^ 

k*l 
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»4 /Jt ■ .2 2 


+ ((B /4 + B u /8)e n +(B , /4)v I -(ITMHBj + 0 n ) 


- <B 3 W /3)($ o + B d (-l) k ) + <B 2 |i 2 /8)B_ T )(l/2) 22 sin 2*. 

k-1 * 


+ ((B y/3)0j I + (B 3 m/6)v i - (B 3 u/6)(Bj + 0^) 


b 

- (B 2 u 2 /4)(f$ o + B d (~l) k ) - B 3 w /6)(l/2) ^ (cos 1 >. - cos 3 

k=l R 


*k> 


.22 


+ ((-B II /4)e o + (B 11/3^ - (Bp/6)Vj I + (B 3 y/6)(e ri - Bj) 

b 

- (B 3 |i/6) 0_) (1/2) 22 (sin 3* - sin *.) + ((B 2 y 2 /8)0 T 

k=l K k 1 

b 

- (B 2 p 2 /8)Bj) (1/2) 22 (cos 2^ - cos 4 + ((-B 2 y 2 /8)e il: 

k s l 


b 

+ (B 2 y 2 /8)0 II )(l/2) 22 (sin 4^ - sin 2^) 


U 

+ K ((u) (.886) (sin 2t)(b)B. + 22 (.2675)p 2 (cos 4*. - 1)0 T 

k=l k 1 

b 

22 (.2675)u sin 4& B TT + (.028)y b -os 2t 8. 

k=l k ii d 

b 

- 22 (.028) (1 - cos*2 . ) (B TT - 8 t ))} 
k=l k II I 
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(The signs of v Q , Vj and Vjj are now changed for positive down 
sign convention.) 

- -(ao/4) { (2B 3 p/ 3) 0 - (B 4 /4 + 3B 2 p 2 /8)0 T 

O I 

- (B 2 1I 2 cos 4t,'8)0, - (B 2 p 2 sin 4t/8)0__ - (B 2 p/2)v 

1 li o 

- (B 3 u/3)B o - (B 4 /4)(8 n - Bl ) - (B 4 /4)v II - (B 2 V 2 /8)B I 
+ (B 2 ji 2 cos 4t/8)B I + (B 2 iT sin 4t/8)B 1I 

- (B 3 w cos 2t/3)B d 4 (B 3 y sin 2t/3)0 d 

+ t ((.886) (y sin 2t)0 d + (.2675)p 2 (cos 4t - ')B T 
+ (.2675)p 2 sin 4t Bjj + (.028)p cos 2t 0 d 

- (.028) (B n - Bj))> (139) 

Similarly 

Cy » -(l/pwQ 2 R 5 )(padl 2 R 4 /2){(B 4 /4 4 B 2 p 2 /4)0 o 

- (B 3 p/3)0 I 4 (B 3 /3)v o 4 (B 3 p/6)v n - (B 3 V i/6)0 I 

b 

- (B 4 /4)(0 o 4 § d (-l) k )) (£ cos * ) 4 ( (2B 3 u /3)0 q 

k=l 

- (B 4 /4 4 3B 2 u 2 /8)$ i 4 (B 2 p/2)v o 4 (B 4 /4) Vii - (B 3 p/3)(0 o 4 8 d (-l) k ) 

- (B 4 /4)(0 n - Bj) - (B 2 p 2 /8)6 I )(l/2) 
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( £ sin 2* fc ) + ( (B 4 /4 + B 2 m 2 /8)0 ii + (B 4 /4)Vj. 
k»l 

- (B 4 /4)(Bj + Bjj) - (B 3 w/ 3)(B 0 + B d (-l) k ) + (B 2 y 2 /8)B n ) 

b 

d/2) ( Y (i + co * 2 \» + C(B 3 i»/3)e II + (B 3 y/6) Vl 

k»l 

- (B^/eKBj + B n ) - (B 2 m 2 /4)(8 o + B d '-l) k ) - B 3 y/6) 

b 

(1/2) ( Y (sln * k + 8in 3 *k» + + (B 3 m/3)0 i 

k«l 


- (B 3 |i/6)v n + (B 3 y/6)(B I1; - Bj) - (B 3 y/6)Bj) (1/2) 

( £ < cos 3 * k + cos * k )) + ((B 2 y 2 /8)0 I - (BV/8)B r ) 
k=l 

b 

(1/2) ( Y (8in 2 * k + a±a A \» + (-(bV/8)0 ii + (B 2 y 2 /8)e tI ) 
k=l 

b b 

(1/2) ( (cos 2tf> k + cos 44» k )) - *( Y f “(*886)11(1 + cos 2 4^^ 

k-1 k=l 

b 

+ Y "(*886) y(l + cos 2* k )(-l) k 8 d + 
k=l 

b b 

Yt "(*886)(ii/2)(3 cos 4» k + cos 3 ^Bj “ (1.070) (y 2 /2) (sin 2iJ» k 

k=l k»l 
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b 

+ sin 4^/2)^ - ^ (*886) (p/2) (sin 3^ + sin <^>8^ 
k-1 

b b 

- £ (1.070) (y 2 /4)(l - cos 4* k )B Ij; + 23 t- 028 )* 1 «i* 2 # k ("D* 8 d 

k-1 k=l 

b 

- £ (.028) (cos 2^ k + 1)(6 I + Bjj))} 

k-1 

The signs of v q , and are changed for positive down sign convention. 

- -(aa/4) ((B 2 u 2 /8)sin 4t e i + (B 4 /4) + B 2 y 2 /8)6 n 

- (B 2 p 2 /8)cos 4t 0 n - (B 4 /4)(l I + B XI + Vj) - (B 3 y/3)& o 

- (B 2 y 2 /8)sin 4t 8 + (B 3 y/3)sin 2t 0 d + (B 3 y/3) 
cos 2t 6 d - (B 2 y 2 /8)B IX + (B 2 y 2 /8)cos 4t 8^ 

+ K 3)yB + (.886)y cos 2t 6, - (.2675)y 2 

O a 

sin 4t 8 - (,2675)y 2 (1 - ccs 4t)8 IT _ - (.028)y 

sin 2t P d - (.028) (Bj. + B XI > )> 


( 140 ) 
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8.3 MULTIBLADE REPRESENTATION OF THE PERTURBATION FLAPPING EQUATIONS, 
INCLUDING REVERSE FLOW 

5 + (Y/2)(CO + CC4 + CN4)§ + (u? + (y/ 2)PS4 sin 4t)S rt 
o ox o 

+ (y/ 2)(-(1/2)CS1 + (CS3/2)cos 4t - (CS5/2)cos 4t 
+ (PC7/2)cos 8t + PC 1/2 + (PC3/2)cos 4t 
+ (PC5/2)cos 4t)8 i + (y/ 2) ((CS3/2)sin 4t + (CS5/2)sin 4t)ftj 
+ (Y/2)((CS3/2)sin 4t + (CS5/2)sin 4t + (PC3/2)sin 4t 

- (PCS /2) sin 4t + (PC7/2)sin 8t)S n + (y/ 2)(CS1/2 - (CS3/2)cos 4t 
+ (CS5/2)cos 4t) Bjj + (y/ 2)(-PS2 sin 2t - PS6 sin 6t)S d 

- ( y/2 ) 'CC2 cos 2t + CC6 cos 6t)S d 

+ (Y/2) (DO + ~C4 cos 4t)v o + (y/ 2)(ES4 sin 4t)Vj (141) 

+ (Y/2)(FO + FC4 cos 4t)v XI 

« (Y/2) (TO + YC4 cos 4t)0j + (y/ 2)(ZS4 sin 4t)0 n 
Bj + 2 S tI - B x + (y/2) (CS3 sin 4t v CS5 sin 4t)6 0 


+ (y/2) (CO + CC2/2 + (CC2/2)cos 4t + CC4 co* 4t 
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+ (CC6/2)cos 4t + (CC6/2)cos StHBj + B n ) + (y/2) ((CC2/2)sin 4t 

- (CC6/2)sin 4t + (CC6/2)sin 8t)(B n - Bj) + (Y/2)(PC1 + 

2 

PC3 cos 4t + PC5 cos 4t + PC7 cos 8t)6 + w, B T 

O XI 

+ (y/2) ((PS2/2)sin 4t + PS4 sin 4t + (PS6/2)sin 4t 
+ (PS6/2)sin + (y/2) (PS2/2 - (PS2/2)cos 4t 

+ (PS6/2)cos 4t - 0’36/2)cos 8t)3 Ii; - (y/2) (CS1 sin 2t 
+ CS3 sin 2t + CS5 sin 6t)B d - (y/2) (PCI cos 2t 

+ PC3 cos 2t + PC5 cos 6t + PC7 cos 608^ 

+ (y/ 2) (DS3 sin 4t + DS5 sin 4t)v o + (y/2)(EC1 + 

EC3 cos 4: + EC5 cos 4t)v^ + (y/ 2)(FS3 sin 4t + FS5 sin 

- (y/2) (YS3 sin 4t + YS5 sin 4t)e r + (y/2)(ZC1 + 

ZC3 cos 4t + ZC5 cos 4t)0jj 

Bjj - 2 Bj - B n + (y/ 2) (CS1 - CS3 cos 4t + CS5 cos 4t.*3 0 
+ (Y/2)((CC2/2)r.in 4t - (CC6/2)sin 4t + (CC6/2)sin 4t 
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+ (CC6/2)sin 8t)(Bj + B n > + (y/2) (CO - CC2/2 - (CC2/2)cos 4t 
+ CC4 cos 4t - (CC6/2)cos 4t - (CC6/2)cos 8t)(B n - Bj) 

+ (y/2) (PC3 sin 4t - PC5 sin 4t + (PC7)sin 8t)B 

o 

+ (y/2) (PS2/2 - (PS2/2)cos 4t + (PS6/2)cos 4t - (PS6/2) 
cos 8t)Bj+to^ Bjj + (y/2)(“(PS2/2)sln 4t + PS4 sin 4t 

- <PS6/2)sin 4t - (PS6/2)sln 8t)B I1; + (y/2)(CSl cos 2t 

- CS3 cos 2t + CS5 cos 6t)3 d - (y/2)(PClsln 2t - PC3 sin 2t 
+ PC5 sin 6t - PC7 sin 6t)B d + (y/2) 

(DS1 - DS3 cos 4t + DS5 cos 4t)v o + (y/2) (EC3 sin 4t 

- EC5 sin 4t>Vj + (y/2)(FSl - FS3 cos 4t + FS5 cos 4t)v u 

- (y/2) (YS1 - YS3 cos 4t + YS5 cos 4t)0. + (y/2)(ZC3 sin 4t 

X 

- ZC5 sin 4t)6 n 

B d - (y/2) (CC2 cos 2t + CC6 cos 6t)B Q - ( y/2) ((CS1/2) 
sin 2t + (CS3/2)sin 2t + (CS5/2)sin 6t) (Bj + B n > 


( 143 ) 
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+ (y/2) ( (CS1/2) cos 2t - (CS3/2)cos 2t + (CS5/2)cos 6t) 

(8n - gj) - (y/2) (PS2 sin 2t + PS6 sin 6t)8 Q 

- (y/2)((PC1/2)cos 2t + (PC3/2)cos 2t + (PC5/2)cos 6t 

+ (PC7/2)cos 6t)8 x - (y/2) ((PCl/2)sin 2t - (PC3/2)sin 2t 
+ (PC5/2)sin 6t - (PC7/2)sin 6t)8 xx + (y/2) (CO + CC4 cos 4t) 
& d + 6 d + (y/2)PS4 sin 4t & d - (y/2)(DC2 cos 2t 
+ DC6 cos 6t)v Q - (y/2) (ES2 sin 2t + ES6 sin 6t)v^ 

- (y/2)(FC2 cos 2t + FC6 cos 6t)v xx 

» - (y/ 2) (YC2 cos 2t)0 x - (y/2)(ZS2 sin 2t + ZS6 sin 6t)0 IX 


The inflow model is same as that obtained without reverse flow (equation 68) 
except the thrust and the moment coefficients C^., and 

C T - (2 ao/b) { (-DTO - DTC4 cos 4t)v Q + (-ETS4 sin 4t) 


Cy are: 


v x + (-FTO - FTC4 cos 4t)v XI + (PTS4 sin 4t + PTS8 

sin 8t)0 + (CTO + CTC4 co- 4t)8 + (-PTS2 sin 2t 

o o 
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- PTS6 sin 6t)3 d + (-CTC2 cos 2t - CTC6 cos 6t)8 d 
+ (1/2) (PTC3 cos 4t + PTO. + CTS1 - CTS3 cos 4t)3j 
+ (1/2) (CTS3 sin 4t)6j + (1/2)(PTC3 sin 4t + 

CTS3 sin 4t)B XI + (l/2)(-CTSl + CTS3 cos 4t)3 n 
+ (YTO + YTC4 cos 4t)0j + (ZTSA Bin 4t)3 XI } 

- - (ao/2b){-(2.DSl + (-2 DS3 + 2 DS5)cos 4t)v Q 

- (2 EC3 - 2 EC5)sin 4t v x - (2 FS1 + (-2 PS3 
+ 2 FS5)cos 4t)v xx “ ((2 PC3 - 2 PC5)sin 4t 

+ 2 PC7 sin 8t)3 Q - ((-2 PCI + 2 PC3)sin 2t + 

2 PC7 sin 6t)8 d + (-PS2 + 2 CO - CC2 + (PS 2 - PS6 
+ 2 CC4 - CC2 - CC6)cos 4t + (PS6 - CC6)cos 8t)3j 
+ ((-2 PS4 + PS2 + PS6 - CC2 + CC6)sin 4t + 

(PS6 - CCS) sin 8t)8 n “ (2 CS1 + (-2 CS3 + 2 CS5) 


(145) 
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cos 4t)0 - ((2 CS1 - 2 CS3)cos 2t + 2 CSS cos 6t)S- 

o a 

+ ((-CC 2 + CC6)sin 4t - CC6 sin 8t)ij + ((-2 CO 
+ CC2) + (-2 CC4 + CC2 + CC6)cos 4t + CC6 cos 8t)8 n 
+ (2 YS1 + (-2 YS3 + 2 YS5)cos 4t)9 ] . + (2 ZC3 - 2 ZC5) 
sin 4t Qjj} 

Cjj - - (aa/2b){-(2 DS3 + 2 DS5)sin 4t v q 

- (2 EC1 + (2 EC3 + 2 EC5)cos 4t)v x - (2 FS3 + 

2 FS5)sin 4t v + (-2 PCI + (-2 PC3 - 2 PC5)cos 4t 

- 2 PC 7 cos 8t)B + ((2 PCI + 2 PC3)cos 2t + 

o 

(2 PC5 + 2 PC7)cos 6t)0 d + ((-2 PS4 - PS2 - PS6 
+ CC2 - CC6)sin 4t + (-PS6 + CC6)sin 8t)B I 
+ ( (-PS2 - 2 CO - CC2) + (PS2 - PS6 - 2 CC4 - CC2 

- CC6)cos 4t + (PS6 - CC6)cos 8t)B I]; + (-2 CS3 - 2 CS5) 
sin 4t B o + ((2 CS1 + 2 CS3)sin 2t + 2 CS5 sin 6t) 8 d 


( 146 ) 
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+ ((-2 CO - CC2) + (-2 CCA - CC2 - CC6)cos 4t - CC6 cos 8t)§ T 
+ ((-CC2 + CC6)sin At - CC6 sin 8t)0 n 

+ (2 YS3 + 2 YS5)sin At 6j + (2 2C1 + (2 ZC3 + 2 ZC5)cos AtJOj-j) 


U47) 



178 


8.4 TABLE OF COEFFICIENTS FOR THE FLAPPING EQUATIONS GIVEN IN 8.3 

KjWO - PTC1 cos \l> + PTC3 cos 3\p + PTS2 sin 2<J> 

+ PTS4 sin 4<J» + PTS6 sin 6* + PTS8 sin 8ip (148) 


y 

PTC1 

PTC3 

PTS2 

PTS4 

PTS6 

PTS8 

.5 

-.2509 

.0157 

-.1 

-.0059 

-.0009 

0 

.6 

-. 3094 

.0271 

-.1379 

-.0103 

-.0016 

.0005 

.7 

-.3723 

.043 

-.1794 

-.0164 

-.0026 

.0009 

.8 

-.4404 

.0642 

-.2234 

-.0245 

-.0039 

.0013 

.9 

-.5146 

.0913 

-.269 

-.035 

-.0056 

.0019 

1.0 

-.5955 

.1252 

-.3151 

-.048 

-.0077 

.002 

C T W 

** CTO + CTC2 cos 2\p + CTC4 cos Aip + CTC6 cos 



+ CTSi sini|i + CTS3 

sin 3i|/ + 

CTS5 sin 5i|/ 

(149) 

u 

CTO 

CTC2 

CTC4 

CTC6 

CTSI CTS3 

CTS5 

.5 

-.3131 

.0106 

-.0015 

-.0002 

.2195 -.005 

.0001 

.6 

-. 319 * 

,0183 

-.0026 

-.0003 

.2551 -.0088 

.0001 

.7 

-.3285 

.0291 

-.0042 

-.0005 

.2863 -.014 

.0001 

.8 

-. 3404 

.0435 

-.0062 

-.0007 

.3122 -.021 

.0002 

.9 

-.3557 

.0619 

-.0089 

-.001 

.3321 -.03 

.0002 

1.0 

-.3748 

O 

QC 

OO 

-.0121 

-.0014 - 

.3454 -.0413 

.0002 
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M-. (9) » YTO + YTC2 cos 2* + YTC4 cos 49 + YTS1 sin 9 

% 

T + YTS3 sin 39 + YTS5 sin 59 + YTS7 sin 79 (150) 


(1 

YTO 

YTC2 

YTC4 

YTS1 

YTS3 

YTS5 

YTS7 

.5 

-.2506 

.2562 

-.0054 

-.4575 

.0489 

.0015 

.0002 

.6 

-.3089 

.3184 

-.0092 

-.5169 

.0667 

.0024 

.0003 

.7 

-.3718 

.3867 

-.0146 

-.5827 

.086 

.0038 

.0004 

.8 

-.4399 

.4619 

-.0217 

-.6535 

.106 

.0056 

.0006 

.9 

-.514 

.545 

-.0308 

QO 

CM 

• 

1 

.1262 

.0079 

.0008 

1.0 

-.5948 

.6373 

-.0421 

-.8048 

.1461 

.0108 

.001 


M n (9) ■ ZTC1 cos 9 + ZTC3 cos 39 + ZTC5 ccs 59 

“ll 

T + ZTS2 sin 29 + ZTS4 sin 49 -I- ZTS6 sin 69 (151) 



ZTC1 

ZTC3 

ZTC5 

ZTS2 

ZTS4 

ZTS6 

.5 

.3578 

-.0516 

-.0017 

.2455 

-.0056 

-.0002 

.6 

.3793 

-.0717 

-.0028 

.3 

-.009 

-.0002 

.7 

.4036 

-.094 

-.0046 

.3576 

-.0148 

-.0003 

.8 

.4304 

-.1182 

-.0068 

.4188 

-.0219 

-.0003 

.9 

.4593 

-.1436 

-.0097 

.4838 

-.031 

-.0003 

1.0 

.49 

-.17 

-.0133 

.5534 

-.042 

-.0003 



^3 
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(♦) » DTO + DTC2 cos 2* + DTC4 cos 4* + DTS1 sin $ 
°T 

+ DTS3 sin 3* + DTS5 sin 5# + DTS7 sin 7* 


V 

DIO 

DTC2 

DTC2 

DTS1 

DTS3 

DTS5 

DTS7 

.5 

.5329 

-.0625 

0 

.3787 

.0209 

.0028 

.0009 

.6 

.5604 

-.09 

0 

.429 

.0302 

.0041 

.0013 

.7 

.5928 

-.1226 

0 

.4708 

.0412 

.0056 

.0018 

.8 

.6303 

-.1601 

0 

.5042 

.0539 

.0074 

.0023 

.9 

.6727 

-.2026 

0 

.529 

.0683 

.0094 

.003 

1.0 

.7201 

-.25 

0 

.545 

.0843 

.0118 

.0037 


M (*) » ETC1 cos + ETC3 cos 3<t + ETC5 cos 5$ + ETC7 cos 7* 

\ 

+ ETS2 sin 2* + ETS4 sin 4* + ETS6 sin 6* 


P 

ETC1 

ETC3 

ETC5 

ETCT 

ETS2 

ETS4 

ETS6 

.5 

.3079 

-.0046 

.0008 

0 

.1123 

.0024 

0 

.6 

.3105 

-.0079 

.0014 

0 

.1320 

.0043 

0 

.7 

.3141 

-.0126 

.0023 

0 

.1503 

.0069 

-.0001 

.8 

.3188 

-.0187 

.0034 

.0004 

.1668 

.0104 

-.0002 

.9 

.325 

-.0266 

.0049 

.0006 

.1813 

.0149 

-.0002 

1.0 

.3326 

-.0365 

.0067 

.0008 

.1935 

.0205 

-.0002 


(152) 


(153) 
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M (<|i) • FTO + FTC1 cos $ + FTC2 cos 2# + FTC4 cos 44 

n* 

+ FTS1 sin 4 > + FTS3 sin 3$ + FTS5 sin 5$ + FTS7 sin 7$ (154) 


V 

no 

FTC1 

FTC2 

FTC4 

FTS1 

FTS3 

FTS5 

FTS7 

.5 

.0551 

.0013 

-.0344 

-.0208 

.4171 

-.0488 

.0053 

.0004 

.6 

.0785 

.0012 

-.0579 

-.0208 

.4171 

-.0488 

.0053 

.0004 

.7 

.102 

.0012 

-.0814 

-.0208 

.4171 

-.0488 

.0053 

.0004 

.8 

.1255 

.0011 

-.105 

-.0208 

.4171 

-.0488 

.0053 

.0004 

.9 

.1490 

.0010 

-.1285 

-.0208 

.4171 

-.0488 

.0053 

.0004 

1.0 

.1726 

.001 

-.1521 

-.0207 

.4171 

-.0487 

.0052 

.0005 
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K(*) • PCI cos ♦ + PC3 cos 3* + PC5 cos 5* + PC7 
+ PS2 sin 2* + PS4 sin 4* + PS6 sin 6# 

|i PCI PC3 PC5 PC7 PS2 

cos 7+ 
PS4 

(155) 

PS6 

.5 

.154 

-.0023 

.0004 

.00004 .0562 

.0012 

-.00005 

.6 

.1863 

-.0048 

.0008 

.00010 .0792 

.0026 

0 

.7 

.2198 

-.0088 

.0016 

.00019 .1052 

.0048 

0 

.8 

.2551 

-.015 

.0027 

.00033 .1334 

.0083 

0 

.9 

.2925 

-.024 

.0044 

.1631 

.0134 

0 





= 0 



1.0 

.3326 

-.0365 

.0067 

.1935 

.0205 

0 

C(*) - 

CO + CC2 cos 2 $ 

+ CC4 cos 

4iJ« + CC6 cos 6$ 



+ CS1 sin 

+ CS3 

sin 3<i + 

CS5 sin 54> 


(156) 


CO 

CC2 

CC4 

CC6 CS1 CS3 

CSS 


.5 

.2233 

-.0026 

.0006 

=0 

.1485 

.0014 

-.0002 

.6 

.2254 

-.0054 

.0013 

0 

.1751 

.003 

-.0004 

.7 

.2288 

-.01 

.0025 

0 

.1993 

.0057 

-.0007 

.8 

.2341 

-.017 

.0043 

0 

.22 

.0097 

-.0012 

.9 

.2418 

-.0273 

.0068 

0 

.2366 

.0157 

-.0019 

1.0 

.2526 

-.0416 

.0103 

0 

.2476 

.024 

-.0028 
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Mg (<J>) - XQ + XC2 cos 2* + XC4 cos 4iji + XS1 sin* + XS3 sin 3* 


V 

XO 

XC2 

XC4 

XS1 

XS3 

.5 

.2781 

-.0562 

-.0006 

.3076 

-.0018 

.6 

.3019 

-.0793 

-.0013 

.3722 

-.0034 

.7 

.3289 

-.1053 

-.0024 

.4393 

-.0061 

.8 

.3589 

-.1336 

-.0042 

.5097 

-.0103 

.9 

.3912 

-.1633 

-.0067 

.5844 

-.0164 

1.0 

.4251 

-.1937 

-.01 

.6647 

-.0248 


%<♦> 

■ Y0 + TC2 cos 2$ -1- YC4 cos 
-1- YS3 sin 3$ + YS5 sin 5$ 

4$ + YS1 sin i> 


W 

Y0 

YC2 

YC4 

YS1 

YS3 

YS5 

.5 

-.1537 

.1547 

-.0008 

-.3061 

.028 

.0004 

.6 

■*“. 186 

.1879 

-.0019 

-.3414 

.029 

.0008 

.7 

-.2195 

.2228 

-.0034 

-.3815 

.0517 

.0013 

.8 

-.2547 

.26 

-.0057 

-.4255 

.0649 

.0022 

.9 

-.292 

.3004 

-.0091 

-.4726 

.0785 

.0035 

1.0 

-.3321 

.3447 

-.0137 

-.5217 

.0919 

.0053 


(157) 


(158) 
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M ($) - ZC1 cos * 4- ZC3 cos 3* + ZC5 cos 4- ZS2 sin 2$ 

II 

+ 2S4 sin 4^ 4- ZS6 sin 6^ (159) 


y 

ZC1 

ZC3 

ZC5 

ZS2 

ZS4 

ZS6 

.5 

.2501 

-.0285 

-.0003 

.153 

-.0009 

0.0 

.6 

.2624 

-.0404 

-.0007 

.1845 

-.0017 

0.0 

.7 

.2765 

-.054 

-.0012 

.2167 

-.0029 

.0002 

.8 

.2923 

-.069 

-.0021 

.25 

-.0048 

.0004 

.9 

.3097 

-.085 

-.0034 

.2842 

-.0075 

.0007 

1.0 

.3284 

-.1022 

-.0051 

.3202 

-.0112 

.0012 


«„ <*) 

= DO + DC2 cos 2\|) + DC4 

cos 4# + DC6 cos 

6* 


o 

4- DS1 sin ip 4- DS2 

sin 3ip 4- DS5 sin 5 $ 


(160) 

M 

DO 

DC2 

DC4 

DC6 

DS1 

DS3 

DS5 

.5 

.313 

-.0106 

.0015 

0 

.2195 

.005 

-.0001 

.6 

.3195 

-.0183 

.0 026 

.0003 

.2551 

.0088 

-.0001 

.7 

.3285 

-.0291 

.0042 

.0005 

.2863 

.014 

-.0002 

.8 

.3404 

-.0435 

.0062 

.0007 

.3122 

.021 

-.0002 

.9 

.3557 

-.0619 

.0089 

.001 

.3321 

.03 

-.0002 

1.0 

.3748 

-.0848 

.0121 

.0014 

.3454 

.0413 

-.0002 
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H. «i) » EC1 cos $ + EC3 cos 34* + ECS cos 5$ 

I 

+ ES2 sin 2$ + ES4 sin 4$ + ES6 sin 6$ (181) 


V 

EC1 

EC3 

EC5 

ES2 

ES4 

ES6 

.5 

.2221 

-.001 

.0003 

.075 

.0005 

-.0001 

.6 

.2228 

-.0021 

.0007 

.089 

.0013 

-.0003 

.7 

.2239 

-.0038 

.0012 

.1025 

.0024 

-.0004 

.8 

.2257 

-.0065 

.0021 

.115 

.0042 

-.00071 

.9 

.2283 

-.0103 

.0034 

.1262 

.0069 

-.0011 

1.0 

.2318 

-.0157 

.0052 

.1359 

.0106 

-.0016 


M u ($) » FO + FC2 cos 2$ + FC4 cos 4$ + FC6 cos 6ifr + FS1 sin 
11 

+ FS3 sin 3tj> + FS5 sin 5$ < l6 2) 


V 

F0 

FC2 

FC4 

FC6 

FS1 

FS3 

FS5 

.5 

.0742 

-.0736 

-.0008 

.0001 

.2245 

-.0018 

.0002 

.6 

.0875 

-.0861 

-.0017 

.0002 

.2280 

-.0036 

.0006 

.7 

.0996 

-.0968 

-.0032 

.0003 

.2337 

-.0064 

.0011 

.8 

.11 

-.1052 

-.0054 

.0005 

.2425 

-.0109 

.002 

.9 

.1182 

-.1105 

-.0088 

.0008 

.2554 

-.0173 

.0033 

1.0 

.1237 

-.1119 

-.0133 

.0012 

.2732 

-.0262 

.0051 
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